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Preface 


This  report  is  the  result  of  an  investigation  into 
the  application  of  linear  stochastic  optimal  estimation 
and  control  techniques  toward  the  solution  of  the  problem 
of  actively  controlling  the  inertial  instrument  test  plat- 
form at  the  Frank  J.  Seiler  Research  Laboratory.  The 
system  is  modeled  as  a linear  system  with  random  (stochastic) 
disturbances.  A forced  separation  concept  is  employed 
in  order  to  investigate  the  effects  of  the  Kalman  filter 
and  the  optimal  controller,  independently.  The  results 
indicate  that  the  optimal  estimation  and  control  system 
is  capable  of  improving  the  performance  of  the  inertial 
instrument  test  platform  but  not  capable  of  meeting  the 
design  specifications  for  the  platform  as  presently  con- 
figured. 
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vided throughout  this  study,  but  also  for  his  guidance 
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Abstract 

The  study  is  directed  toward  the  analysis  and  im- 
plentation  of  an  ootimal  estimator  (Kalman  filter)  and 
an  optimal  regulator  to  provide  active  control  of  the 
inertial  instrument  test  platform  at  the  Frank  J.  Seiler 
Research  Laboratory.  The  design  specifications  are  to 
maintain  angular  position  within  ±1.0x10  ^ arsceconds 
and  angular  rate  with  ±1.667x10*^  arcseconds/seconds . 

A forced  separation  concept  is  utilized  to  allow 
the  independent  evaluation  of  the  Kalman  filter  and  the 
optimal  regulator.  Optimal  and  suboptimal  Kalman  filter 
models  are  developed  and  evaluated  at  physically  realizable 
sampling  rates.  A general  optimal  estimation  and  control 
algorithm  is  developed  and  a proposed  sequence  of 
algorithm  computations  is  presented. 


ANALYSIS  AND  IMPLEMENTATION  OF  OPTIMAL 
ESTIMATION  AND  CONTROL  FOR  THE  FJSRL 
SEISMIC  ISOLATION  Pi^ATFORM 

I . Introduction 

The  Seismic  Isolation  Platform  at  the  United  States 
Air  Force  Academy's  Frank  J.  Seiler  Research  Laboratory 
(FJSRL)  is  an  inertial  instrument  test  platform  designed 
to  provide  a high  degree  of  isolation  from  environmental 
disturbances.  An  active  control  system  is  needed  +o  pro- 
vide the  isolation  required  for  testing  and  evaluating 
highly  advanced  inertial  components  and  systems  (Ref  1) . 

Background 

As  depicted  in  Figures  1 and  2,  the  inertial  instru- 
ment test  platform,  as  viewed  from  above,  is  a 25  feet  by 
25  feet  square  with  nine  circular  test  tables  extending 
approximately  2 feet  from  the  top  surface  of  the  platform. 
Viewed  from  below,  the  platform  is  cruciform  shaped.  The 
platform  is  constructed  of  steel  reinforced  concrete,  is 
9 feet  high,  and  weighs  approximately  450,000  pounds.  The 
platform  is  located  beneath  a false  floor  in  the  laboratory 
and  the  test  tables  protrude  through  holes  in  this  floor. 
The  platform  is  supported  by  twenty  pneumatic  cylinders 
that  essentially  float  the  platform  a fraction  of  an  inch 


above  the  base  slab.  The  pneumatic  cyclinders  rest  on 
red  oak  blocks  that,  in  turn,  rest  on  the  base  slab. 

The  base  slab  rests  on  a compacted  aggregate  fill  base 
that  is  designed  to  minimize  the  coupling  of  vibrations 
in  order  to  separate  the  platform  base  from  the  building 
foundation.  The  pneumatic  cylinders  are  arranged  as  de- 
picted in  Figure  3.  The  inner  twelve  cylinders  regulate 
the  height  of  the  platform  (referenced  to  the  base  slab) 
and  the  outer  eight  cylinders  are  used  in  a push-pull  con- 
figuration to  regulate  angular  motion  about  the  horizontal 
axis. 

The  system  consisting  of  the  platform  and  the  pneu- 
matic cylinders  has  a natural  frequency,  as  measured  by 
FJSRL  in  January  1977,  of  1.3Hz  and  acts,  effectively,  as 
a passive  isolation  system  (low  pass  filter)  for  distur- 
bances above  1.3Hz  (Ref  3).  Unfortunately,  many  of  the 
disturbances  of  interest,  e.g.  earthquakes,  ocean  waves, 
and  barometric  pressure  variations,  have  frequencies  below 
1.3Hz  and,  therefore,  an  active  control  system  is  required 
to  isolate  the  platform  from  these  disturbances. 

In  order  to  provide  active  control  of  the  platform, 
a combination  of  tiltmeters,  angular  motion  sensors,  and 
actuators  are  attached  to  the  platform  as  depicted  in 
Figure  4.  The  eight  angular  motion  control  cylinders 
(Fig  3),  in  combination  with  the  tiltmeters  (on  the  sur- 
face of  the  platform) , provide  closed-loop  control  of 
angular  motion  about  the  platform's  center  cf  gravity. 
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0 Angular  Motion  Control  Cylinders 
o Height  Control  Cylinders 


Figure  3.  Location  of  Pneumatic  Cylinders  (Ref  2\  17) 

I | 

Additional  damping  control  is  provided  by  a second  control 
loop  consisting  of  angular  motion  sensors  (seismometers) 
and  electromagnetic  one-dimensional  dampers  (shakers) . 

The  shakers  are  attached  at  the  four  corners  of  the  plat- 
form at  the  aDproximate  level  of  the  center  of  gravity. 

El 

A separate  height  control  system  consists  of  the 
\ i 

twelve  inner  pneumatic  cylinders,  and  sensors  that  measure 

. 

•j  the  distance  between  the  base  slab  and  the  bottom  of  the 

platfcrm.  Since  the  height  control  is  not  significantly 

it 

affected  by  external  disturbances  (Ref  4:3),  it  is  not 
mentioned  further  in  this  study. 

The  platform  was  constructed  with  the  cruciform  bottom 


in  order  to  locate  the  center  of  gravity  at  the  level  of 
the  pneumatic  cylinders.  In  this  way,  the  coupling  be- 
tween the  various  modes  of  motion  are  minimized. 
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From  the  beginning,  the  problem  of  interest  has  been 
that  of  designing  and  implementing  controllers  to  provide 
the  desired  degree  of  isolation  (angular  position  within 
±1.0  x 10-3  arcsecond  and  angular  rate  within  ±1.667  x 
10~5  arcseconds/second) . Early  investigations  and  attempts 
at  stabilizing  the  platform  were  primarily  analog  control 
systems  (Ref  5-8) . The  latest  analog  design,  although 
successful  in  meeting  the  specification  for  angular  posi- 
tion, was  unable  to  meet  the  angular  rate  requirements 
(Ref  8) . 

More  recently,  attempts  have  been  made  to  provide 
digital  control  to  the  platform  (Ref  1) , including  attempts 
at  implementing  various  Finite  Impulse  Response  (FIR)  and 
Infinite  Impulse  Response  (HR)  filters  at  FJSRL  (Ref  9)  . 
Digital  control  is  desirable  because  it  is  less  affected 
by  noise  and  other  disturbances  than  analog  systems,  per- 
mits the  use  of  sensitive  control  elements  with  relatively 
low  energy  signals , and  has  phase  characteristics  that 
cannot  be  duplicated  by  an  analog  system  (Ref  9) . Digital 
control  is  more  flexible  than  analog  control  because  changes 
in  sensors  or  actuators  can  often  be  incorporated  with 
only  software  modifications.  An  additional  benefit  of 
digital  controllers  is  their  adaptability  to  the  solution 
of  stochastic  optimal  estimation  and  control.  It  is  this 
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latter  trait  of  digital  control  that  is  most  useful  for 
this  investigation. 

To  date,  the  attempts  at  digital  control  of  the  plat- 
form have  not  been  successful.  Although  some  design  work 
has  demonstrated  the  feasibility  of  digital  control  (Ref  1: 

9) , the  actual  implementations  of  these  designs  have  not 
been  successful.  The  failure  of  the  implemented  controllers/ 
filters  to  meet  the  theoretical  performance  levels  is  due 
primarily  to  the  effects  of  slow  sampling  rates,  finite 
wordlengths , the  uncertainties  in  the  system  models,  pro- 
cess noise,  and  the  sensitivities  of  the  sensors  employed. 

In  the  previous  digital  attempts  (Ref  1:  9),  the  platform 
was  modeled  as  a completely  deterministic  system. 

In  1976,  two  investigations  were  completed  in  which 
the  platform  was  modeled  as  a stochastic  system  (Ref  2:  4). 

Stochastic  modeling  of  the  system  permits  the  application 
of  optimal  estimation  and  stochastic  control  methods  to 
the  problem  of  isolation  of  the  test  platform.  Optimal 
estimation  is  desirable  because  it  allows  consideration 
of  the  stochastic  (random)  characteristics  of  the  platform 
system,  including  system  modeling  errors,  the  randomness 
of  the  environmental  disturbances,  the  process  noise,  and 

t 

the  error  sources  attributed  to  the  sensors. 

n 

Statement  of  the  Problem 

The  purpose  of  this  study  is  to  analyze  and  implement 
an  optimal  estimator  and  controller  for  use  at  FJSRL  for 

i 

I 

j 

■ I 
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the  active  control  of  the  seismic  isolation  platform. 

The  specifications  are  dictated  by  the  performance 
criteria  for  new  generation  inertial  components  (Ref  5) . 
These  criteria  require  that  the  angular  position  (tilt) 
of  the  platform  be  maintained  within  ±0.001  arcseconds 
and  the  angular  rate  (velocity)  be  maintained  within 
±0 . 001  arcseconds  per  minute  (1.667xl0~5  arcseconds  per 
second) . In  addition,  the  controller  must  be  effective 
in  the  frequency  band  of  0-20  Hz  (Ref  8)  with  a step  in- 
put (for  testing  purposes)  of  2.5  foot-pounds.  A step 
input,  applied  directly  to  the  top  surface  of  the  block, 
is  used  for  testing  because  it  is  easily  modeled  in  the 
s-domain  and  also  in  the  discrete  (z)  domain.  The  2.5 
foot-pound  step  input  was  chosen  because  it  has  been  demon- 
strated that  it  approximates  the  disturbance  caused  by 
moderate  environmental  disturbances  on  the  passive  plat- 
form (Ref  2:20). 

Assumptions 

The  system  is  considered  to  be  linear  over  the  fre- 
quency band  of  0-20  Hz  and  it  is  assumed  that  the  construc- 
tion of  the  platform  and  the  locations  of  the  sensors  is 
such  that  the  coupling  of  the  modes  of  motion  is  minimized 
to  the  extent  that  they  can  be  considered  separately.  The 
linearity  assumption  has  been  supported  by  previous  studies 
(Ref  8)  and  is  justified  by  the  fact  that  only  small  pertur- 
bations occur  when  the  platform  is  being  controlled. 

8 


Structural  resonances  above  20  Hz  (the  total  number 


is  unknown)  are  disregarded.  It  is  assumed  that  these 
structural  modes  will  not  have  a significant  impact  on 
the  performance  of  the  system  due  to  natural  passive  damping 
(Ref  8).  If  it  is  determined  that  they  do  have  an  effect, 
the  model  can  be  redesigned  to  take  them  into  account. 
Pseudo-noise  can  be  added  to  the  system  model  to  represent 
the  inaccuracies  in  the  system  model  resulting  from  struc- 
tural resonances. 

The  process  and  measurement  noises  can  be  modeled 
as  White  and  Gaussian.  The  White  assumption  is  valid  due 
to  the  wideband  characteristics  of  the  noise  and  the  very 
narrow  band  characteristics  of  the  system  (wideband  noise 
driving  a limited  bandwidth  system) . The  Gaussian  assump- 
tion is  based  on  the  Central  Limit  Theorem  since  the  pro- 
cess and  measurement  noises  are  composed  of  several  inde- 
pendent additive  noises  (Ref  10) . The  process  noise  and 
measurement  noise  are  considered  to  be  independent  (Ref  4:6). 
This  assumption  is  based  on  the  fact  that  the  sensors  are 
separated  and  the  measurement  process  does  not  corrupt  the 
state  being  measured. 

A PDP-11/03  minicomputer  has  been  designated,  by  FJSRL , 
for  use  in  implementing  the  digital  control  system  for  the 
seismic  isolation  platform.  The  central  computing  system  in 
the  PDP-11/03  is  the  LSI-11  minicomputer  board.  The  filter/ 

A 

controller  algorithm  developed  in  this  study  is  implemented 
in  LSI-11  compatible  assembly  code. 
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Design  Approach 

The  approach  is  defined  as  the  solution  of  a linear 
stochastic  estimation  and  control  problem.  A linear 
stochastic  system  is  described  by  the  matrix  equations 


Fx ( t)  + Lu ( t)  + Gw ( t) 


z(tk) 


Hx(tk)  + v(tk) 


where  F,  L,  G,  and  H are  time  invariant  matrices  derived 
from  the  system  transfer  functions.  The  column  vector 
x(t)  is  the  state  variable  vector,  u(t)  is  the  control 
input  vector,  w(t)  represents  the  environmental  disturbances 
to  the  noise  shaping  filter,  £(tk)  is  the  sampled  measure- 
ment vector,  and  v(tk)  represents  the  sampled  measurement 
noise. 

The  objective  of  the  linear  stochastic  control  problem 
is  to  find  the  discrete- time  control  input 


u(tk)  = -C(tk)£(tk) 


that  minimizes  the  quadratic  performance  index 


J * E{l/2xT(tN+1)Vfx(tN+1)  + 


I l/2[xT(tk)V(tk)x(tk)  + uT(tk)U(tk)u(tk)  ]}  (4) 


where  is  a discrete-time  estimate  of  the  state  vector, 

x(tk) , is  the  discrete  control  input  vector  (put 

through  a zero-order  hold,  ZOH) . E is  the  expected  value 
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operator,  and  tN+^  is  the  terminal  time.  The  matrices  Vf, 

V,  and  U are  weighting  matrices  that  assign  performance 
costs,  over  the  entire  time  interval  of  interest,  due  to 
the  terminal  states,  the  present  states,  and  the  control 
inputs,  respectively.  The  performance  of  the  optimal  con- 
troller is  a function  of  the  values  chosen  for  these 
weighting  matrices.  Only  through  the  proper  selection  of 
these  weighting  matrices  can  u(t^)  meet  the  design  speci- 
fications. 

With  the  assumptions  previously  stated,  employment  of  the 

Separation  Principle  was  considered  to  simplify  the  design 

and  analysis  of  the  system.  The  Separation  Principle  states: 

The  optimal  stochastic  controller  for  a linear  system 
driven  by  white  Gaussian  noise,  subject  to  a quadratic 
cost  criterion,  consists  of  an  optimal  linear  Kalman 
filter  cascaded  with  the  optimal  feedback  grain  matrix 
of  the  corresponding  deterministic  optimal  control 
problem  (Ref  11:11-16). 

The  Separation  Principle  is  depicted  in  Figure  5. 

In  the  investigation,  a "forced  separation"  of  the 
Kalman  filter  and  the  optimal  controller  is  invoked,  based 
on  engineering  judgement,  to  take  advantage  of  some  of  the 
physical  properties  of  the  system  elements.  Because  of 
the  inherent  dynamics  of  the  system  elements  (the  pneumatic 
actuator  has  a settling  time  of  20  seconds  and  the  electro- 
magnetic actuator  has  a natural  frequency  of  26  Hz)  and, 
because  the  pneumatic  actuator  is  actually  part  of  the  plat- 
form support  system,  the  states  associated  with  the  pneu- 
matic actuator  are  not  included  in  the  system  state  space 
models  used  in  the  analyses  of  the  Kalman  filter.  In 
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addition,  the  states  associated  with  the  sensors  are  not 
included  in  the  system  model  used  in  the  development  of 
the  optimal  controller  because  these  states  contain 
stochastic  properties  (measurement  noises)  and  it  was  con- 
sidered desirable  to  eliminate  stochastic  effects  from 
the  development  and  analysis  of  the  optimal  controller. 

The  concept  of  a "forced  separation"  is  necessary  because, 
although  the  system  being  investigated  meets  the  White 
Gaussian  noise  and  quadratic  cost  function  criteria,  the 
Separation  Principle  does  not  apply,  in  the  strict  sense, 
because,  in  this  investigation,  the  states  represented  by 
the  Kalman  filter  model  are  not  the  same  states  as  those 
represented  by  the  optimal  controller  model.  By  utilizing 
the  concept  of  "forced  separation",  the  optimal  estimation 
and  control  problem  is  separated  into  the  design  and  analysis 
of  a Kalman  filter  independent  of  the  design  and  analysis 
of  the  corresponding  optimal  controller.  Without  the 
"forced  separation",  the  development  and  analysis  of  the 
optimal  estimation  and  control  problem  would  be  far  more 
complex,  if  at  all  tractable. 

The  designs  of  a Kalman  filter  and  a deterministic 
optimal  controller  were  investigated  by  Richard  Brunson 
and  Martin  J.  Burkhart,  respectively  (Ref  2:  4).  The 
approach  taken  n cnis  investigation  is  based,  in  part, 
on  their  efforts.  Both  investigations  were  based  on  a 
sampling  rate  of  200  Hz  (based  on  engineering  judgment, 
using  five  times  the  Nyquist  frequency*)  and  involved, 
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The  investigation  includes  the  effects  of  conversion  quan- 
tization (A/D  and  D/A),  coefficient  quantization,  and 
roundoff,  truncation,  and  overflow  due  to  arithmetic  oper- 
ations. In  addition  the  numerical  precision  problems 
associated  with  the  Kalman  filter  are  discussed,  along 
with  the  appropriate  techniques  for  decreasing  the  effects 
of  these  problems. 

Previously  derived  component  transfer  functions  (Ref 
2:  4)  are  used  to  develop  a system  "truth  model",  in  state 
variable  form,  as  the  basis  for  designing  an  optimal  Kalman 
filter.  The  optimal  filter  is  analyzed  for  the  effects  of 
sampling  rate  and  finite  wordlength.  In  addition,  a sensi- 
tivity analysis  is  performed  (using  covariance  analysis 
techniques)  to  demonstrate  the  effects  of  the  various  noise 
sources  separately.  Using  various  simplifying  assumptions, 
four  suboptimal  Kalman  filter  models  (reduced  order  filters) 
are  developed  for  comparison  and  possible  implementation. 

Using  component  transfer  functions,  an  optimal  con- 
troller  is  developed  under  the  assumption  that  it  is  re 
ceiving  "perfect"  information  from  the  Kalman  filter,  i.e. 
exact  knowledge  of  the  entire  state.  The  design  of  the 

optimal  controller  is  based,  primarily  on  the  linear  quad- 
i 

ratic  full-state  feedback  controller  developed  by  Burkhart 
(Ref  2) . This  model  is  derived  for  a sampling  rate  of 
200  Hz. 

A general  estimation  and  control  algorithm  is  designed 
and  each  filter  model  investigated  (in  cascade  with  the 
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appropriate  optimal  controller)  to  determine  which  filter 
provides  the  best  results.  This  is  accomplished  by  deter- 
mining the  computation  time  involved  for  each  combination 
and,  thus,  the  maximum  sampling  rate  possible  for  each 
combination.  Each  Kalman  filter  model  is  then  tuned  at 
this  maximum  sampling  rate  and  the  results  compared  by 
utilizing  a covariance  analyses.  A brief  description  of 
considerations  for  implementing  the  algorithm  is  presented 
as  a baseline  for  future  investigations. 

Organization  of  Report 

The  report  is  divided  into  five  chapters.  In  Chapter 
II,  the  "truth  model"  Kalman  filter  is  derived  and  analyzed 
and  various  suboptimal  filters  are  developed.  Prominent 
in  the  analysis  are  the  effects  of  sampling  rate,  finite 
wordlength,  and  noise  sources. 

In  Chapter  III,  the  optimal  controller  is  discussed 
in  terms  of  the  expected  results  and  previous  conclusions. 

In  Chapter  IV,  a general  optimal  estimation  and  con- 
trol algorithm  is  developed.  Each  Kalman  filter  is,  in 
turn,  combined  with  the  appropriate  optimal  controller 
and  analyzed.  Prominent  in  the  analysis  are  the  filter 
tuning  process  and  the  effects  of  sampling  rate.  The 
results  of  this  analysis  are  compared  and  a "best  cut" 
combination  selected.  The  effects  of  quantization,  con- 
version, computational  delay,  and  specific  characteristics 
of  the  minicomputer  (LSI-11)  are  discussed. 
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II.  Kalman  Filter  Development 


The  development  of  the  optimal  Kalman  filter  is  based 
on  the  state  variable  (truth)  model  of  the  seismic  plat- 
form system  composed  of  the  platform  dynamics,  the  sen- 
sors, and  the  process  and  measurement  noise  sources.  In 
this  chapter,  the  system  truth  model  is  developed  and  the 
optimal  Kalman  filter  is  designed  and  analyzed.  The 
Kalman  filter  based  on  the  system  "truth"  model,  i.e.  the 
optimal  Kalman  filter,  is  used  as  a "benchmark"  to  deter- 
mine the  best  performance  that  can  be  expected  from  the 
optimal  estimation  of  the  state  of  the  seismic  isolation 
platform  system.  The  effects  of  finite  wordlength,  sampling 
rate  and  noise  levels  are  investigated.  In  addition,  by 
making  various  simplifying  assumptions,  four  sub-optimal 
Kalman  filters  are  developed  based  on  reduced-order  system 
models.  In  Chapter  IV,  these  suboptimal  filters,  and  the 
appropriate  optimal  filter,  are  combined,  in  turn,  with 
the  optimal  controller  developed  in  Chapter  III.  Each  fil- 
ter/controller combination  is  analyzed  and  compared  for 
implementation  in  the  form  of  a Linear  Quadratic  Gaussian 
(LQG)  controller. 

Justification  for  Kalman  Filter 

As  mentioned  previously,  attempts  at  controlling  the 
platform  using  deterministic  digital  control  methods  failed 
to  provide  the  specified  control  of  the  seismic  isolation 
platform.  Deterministic  approaches  do  not  include  con- 
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siderations  of  the  effects  of  uncertainties  in  the  system 
model,  random  effects  of  process  and  measurement  noises, 
or  the  sensitivities  of  the  sensors  employed. 

Although  the  best  model  available  for  the  seismic 
isolation  platform  is  used  for  this  study,  there  remain 
uncertainties  that  must  be  considered.  There  is  some 
question  as  to  the  homogeneity  of  the  platform  structure 
and  this  leads  to  uncertainties  about  the  exact  center  of 
gravity.  The  assumption  of  the  decoupling  of  the  modes 
of  motion  might  not  be  valid.  Uncertainties  about  the 
bending  modes  and  resonant  frequencies  add  to  the  overall 
inaccuracy  of  the  system  "truth  model".  In  addition,  there 
is  process  noise,  measurement  noise,  and  possible  biases 
in  the  measurement  sensors  that  must  be  considered.  A 
difficulty  in  determining  the  current  state  of  the  plat- 
form angular  rate  arises  because  there  is  no  direct  measure- 
ment of  this  state;  rather,  it  is  determined,  indirectly, 
from  the  measurements  of  the  platform  tilt  and  rate.  Also, 
the  sensors  are  most  accurate  in  different  frequency  bands. 
The  tiltmeter  is  effective  in  the  0-1  Hz  range,  and  the 
angular  rate  sensor  is  more  accurate  in  the  higher,  1-20 
Hz  range.  All  of  these  factors  contribute  to  inaccuracies 
in  the  determination  of  current  system  states.  A stochastic 
estimator  is  required  to  filter  the  effects  of  the  "noise” 
due  to  the  above  factor.  Under  the  assumption  of  linearity 
and  white  Gaussian  noise  sources,  it  can  be  shown  that  a 
Kalman  filter  is  the  "optimal"  estimator  of  the  current 


18 


state  of  the  system  (Ref  11:1-44). 


A Kalman  filter  is  a recursive  data  processing  algo- 
rithm that  uses  information  about  the  system  dynamics, 
initial  states,  and  statistics  of  the  process  and  measure- 
ment noise  to  generate  an  "optimal"  estimate  of  the  current 
states  of  the  system.  A Kalman  filter  attempts  to  minimize 
the  uncertainties  and  random  noise  from  a system  through 
the  application  of  all  current  measurements  combined  with 
the  assumed  statistics  of  the  system.  Since  it  is  recur- 
sive, it  is  not  required  that  all  past  information  be 
remembered  (stored  in  the  computer) ; the  estimation  and  co- 
variance  from  the  previous  update  are  sufficient  statistics. 

The  optimal  estimation  problem  is  separated  from  the 
optimal  control  problem,  in  this  investigation,  by  invoking 
the  forced  separation  concept  described  in  Chapter  I. 

This  permits  the  independent  design  and  analysis  of  the 
optimal  estimator  (Kalman  filter)  without  consideration  of 
the  optimal  control  gain  (addressed  in  Chapter  III) . The 
outputs  of  the  Kalman  filter  (state  estimates)  become  the 
current  state  inputs  to  the  optimal  controller  when  the 
controller  and  estimator  are  recombined  (Chapter  IV) . 

Brunson  investigated  the  feasibility  of  employing  a 
Kalman  filter  to  improve  the  estimates  of  the  system  states 
of  the  seismic  platform  (Ref  4) . Because  the  design  was 
based  on  early,  less  accurate,  system  models  and,  because 
the  implementation  was  analyzed  for  a HP-21MX  minicomputer, 
the  results  are  used  in  this  investigation,  primarily,  as 
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background  reference  and  as  a basis  for  possible  design 
techniques. 


Development  of  System  Truth  Model 

The  seismic  isolation  platform  system  used  for  the 
design  of  the  Kalman  filter  consists  of  the  platform, 
sensors  (tiltmeters  and  seismometers) , and  noise  sources. 
The  platform  system  can  be  modeled  as  a linear  stochastic 
system  and  can  be  represented,  in  state  variable  form, 
by  the  equation 

*(t)  = F ( t)  x ( t)  + G (t)  w (t)  + L (t)  u (t)  (5) 

where 

F (t)  system  matrix 

G ( t)  matrix  of  states  corrupted  by  process  noise 
L ( t ) control  matrix 
x(t)  state  vector 
w(t)  process  noise  vector 
u(t)  control  input  vector 
The  measurement  equation  for  the  system  is  given  by 

^(t^)  = H(tk)x(tk)  + v(tk)  (6) 

where 

^(t^)  measurement  vector 

H ( t,  ) measurement  matrix 
k 

x(tk)  state  vector 

v(t^)  measurement  noise  vector 
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The  matrices  F(t),  G(t),  L(t),  and  H(t)  are  assumed 
to  be  time-invariant  and  are,  therefore,  expressed  as 
F,  G,  L,  and  H through  the  remainder  of  this  thesis. 

The  truth  model  is  depicted,  in  block  diagram  form,  in 
Figure  6. 

The  process  noise,  w(t) , and  the  measurement  noise, 
v(t^),  are  considerad  to  be  independent  and,  therefore, 
uncorrelated. 

The  process  noise  is  described  as  zero-mean,  White 
Gaussian  noise.  The  process  noise  covariance  matrix,  Q(t), 
is  given  by 


E [w(t) w (t  + t) ] = Q6t 


(7) 


The  measurement  noise  is  described  as  zero-mean,  white 
Gaussian  noise  and  the  measurement  noise  from  the  tilt- 
meter  is  assumed  to  be  independent  of  the  measurement 
noise  from  the  seismometer.  The  sampled  measurement  noise 
covariance  matrix,  R(t^),  is  given  by 


E[v(ti)v(tk)T] 


R(t.) 


i-k 

i*k 


(8) 


In  the  following  sections,  the  F,  G,  L,  and  H matrices 
are  derived  from  the  system  transfer  functions  (s-domain) . 

System  Dynamics  Model  (Platform) . The  platform  trans- 
fer function,  for  a torque  input  and  horizontal  angular 
output,  is 


21 


A 


Figure  6.  System  Truth  Model 
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Vb 
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Arcseconds 
Foot-pounds  (9) 


where 

5b  = 0.05 

a)b  = 7 rad/sec 

Kb  = 0.044  arcsecond/f t-lb  (Ref  2:16) 

Solving  the  equation,  for  the  given  variables, 

6 _ 2.156  Arcseconds 

_ = 

T s + 0.07s  + 49  Foot-pounds 


yields 

(10) 


where  0 is  the  angle  of  the  platform  referenced  to  local 

level  and  T is  the  external  torque  applied  to  the  platform. 

T is  composed  of  torques  from  the  process  noise,  represented 

by  T , pneumatic  actuator,  u , and  electromagnetic  actuator 
w * p 

(shaker),  ug.  Therefore, 


T = T + u + u„  (Foot-pounds) 
w p s 


(ID 
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The  transfer  function  is  converted  to  state  variable 


form  using  phase  variable  techniques  because  the  states 
of  interest,  i.e.  angular  position  (0)  and  angular  rate 
(0)  are  represented  directly  in  the  resulting  form.  Letting 


0 

0 


(Arcseconds) 
(Arcseconds/Second ) 


(12) 

(13) 


and  cross-multiplying  the  transfer  function,  the  state 
equation  become 


x-  = -4 9x.  - 0. 7x_  + 2.156T  + 2.156u  + 

2 12  W p 

+ 2 . 156u_ 


Tiltmeter.  The  tiltmeter  transfer  function,  from 
Reference  4:12,  is 


(14) 


(15) 


Vx  19.84(1  - 0. 1531s  ) 


Volts 


~2 


■(16) 


s + 12.56s  + 77.6s  + 198.4  Arcsecond 


where  is  the  tiltmeter  output.  Based  on  results  from 
Brunson  (Ref  4:48),  the  tiltmeter  model,  and  subsequent 
component  models  are  derived  in  physical  variable  form. 
Therefore,  letting 


(17) 


the  state  equations  become 


x^  = -3.037504x^  - 12.56x2  + 


[ 18) 


-77.6x2  + x^ 


(19) 


* 19.84x^  - 198.4x2 


and  the  measurement  equation  is 


(20) 


(21) 


Angular  Motion  Sensor.  The  angular  motion  sensor 
(seismometer)  transfer  function  is  (Ref  4:12) 


W_2 

s^  Volts 

6 

2 

s + 2s  + 1 Arcsecond 

(22) 

where  V2  is  the 

state  equations 

seismometer  output.  The  physical 

are 

variable 

X6  = 

2xl  " 2x6  + x7 

(23) 

x?  = xL  - x6 

The  measurement  equation  is 

(24) 

Z2  = 

*1  - x6  * V2 

(25) 

Process  Noise  Filter.  The  process  noise  represents 
disturbance  torques  from  external  sources,  such  as  earth 
seismic  activity,  that  are  input  to  the  platform.  The 
process  noise  is  modeled  as  a time  correlated  Gaussian 
noise  plus  a white  noise  (Ref  10) . A noise  shaping  filter 
is  required  in  order  to  provide  noise  with  the  desired 
power  spectral  density  properties  through  the  range  of 
0-20  Hz.  Brunson  developed  a third-order  approximation 
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for  the  noise  shaping  filter,  called  the  process  noise 
filter  in  this  study  (Ref  4:14).  The  transfer  function 
is 


C 

1 

(s  + 125) 3 

(26) 

which 

yields  the 

: state  equations 

x • 

00 

II 

-125Xg  + xg 

(27) 

x9  = 

“125Xg  + x10 

(28) 

x10  = 

_125x10  + W1 

(29) 

where 

X8 

is  the 

output  from  the  process 

noise  filter  that 

is  input 

to  the 

platform  in  the  form  of 

disturbance  torque, 

called  Tw  (see  platform  model  description),  and  w^  represents 
the  white  noise  input  to  the  filter. 

Pneumatic  Actuator  Dynamics . The  pneumatic  actuator 
dynamics  are  associated  with  the  control  of  the  platform 
and  are  omitted  during  analysis  of  the  Kalman  filters 
(forced  separation). 

Shaker  Dynamics . In  the  frequency  range  of  interest 
(0-20  Hz),  the  shaker  has  essentially  no  dynamics.  Figure 
7 represents  the  frequency  response  of  the  shaker  actuator. 
For  purposes  of  analysis,  the  shaker  transfer  function  is 
replaced  by  a constant  gain. 

Combined  State  Equations . The  full  set  of  state  equa- 
tions representing  the  truth  model  is 
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X1 

" 

X2 

X2 

= 

-49x1  - 0.7x2  + 2 . 156Xg 

X3 

= 

-3.037504x1  - 12.56x3  + 

X4 

-77.6x3  + x^ 

X5 

= 

19.84x1  - 198. 4x3 

X6 

3 
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x7 

= 
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X • 
00 

= 
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x9 

= 
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3 
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The  F matrix  is 

0 10 
-49  -0.7  0 

-3.037504  0 -12.56 


0 

19.84 

2 

1 

0 

0 
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0 -77.6 

0 -198.4 

0 0 
0 
0 
0 
0 


0 

0 

0 

0 


0 

0 

1 

0 

0 

0 

0 

0 

0 

0 


0 

0 

0 

1 

0 

0 

0 

0 

0 
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0 
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and  the  G and  L matrices  are 


(30) 


K31) 
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The  control  input  is 


The  control  input  u^,  associated  with  the  pneumatic 
actuators,  is  generated  by  a deterministic  position  feed- 
back controller  that  is  used  to  directly  counter  any 
torque  imbalance  resulting  from  unsymmetrical  positioning 
of  items  on  the  platform.  The  control  input  ug,  associated 
with  the  electromagnetic  shaker  activators,  is  generated 
by  an  optimal  state  feedback  controller.  The  optimal  con- 
troller is  used  to  regulate  the  position  feedback  controller 
and  also  to  control  the  angular  rate  of  the  platform.  The 
Kalman  filter  receives  input  from  the  sensors  and  the 
optimal  controller  (Fig  8)  and  the  control  input  u^  is 
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Figure  8.  Filter/Controller  Configuration 

not  considered  in  the  design  and  analysis  of  the  Kalman 
filter.  This  is  addressed  further  in  Chapter  III  in  the 
discussion  of  the  design  of  the  controllers.  Because  of 
the  "forced  separation",  ug  is  the  only  control  input  of  in- 
terest in  the  Kalman  filter  analysis,  the  Lu(t)  term  becomes: 

0 

2.156 
0 
0 
0 
0 
0 
0 
0 
0 
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Finally,  the  system  measurement  matrix 


is 


Kalman  Filter  Mode  1 

A Kalman  filter  is  a predictor-corrector  type  of 
estimator  that  uses  a conditional  probability  density, 
conditioned  on  the  actual  measurements,  to  describe  the 
probabilities  of  possible  system  states.  The  conditional 
probability  density  is  a function  of  the  system  dynamics, 
initial  states,  and  the  assumed  statistics  of  the  distur- 
bance noises.  Since  the  conditional  probability  density 
function  itself  is  Gaussian,  it  is  completely  described 
by  the  first  and  second  order  statistics,  i.e.  conditional 
mean  and  covariance.  The  general  Kalman  filter  equations 
that  represent  these  statistics  are  divided  into  two 
functions;  those  that  propogate  (predict)  the  conditional 
mean  (optimal  estimate)  and  the  covariance,  and  those 
that  update  (correct)  the  optimal  estimate  and  covariance 
at  measurement  sample  times. 

The  propagation  equations  (in  discrete  form)  are 


*(tk} 


$ (tk,tk-l)  -(tk-l] 


+r (tk,tk-l)-(tk-l) 


(36) 


and 

p(t-) 


+ 


/k  |f’(tk,T)G(T)Q(T)GT(T)  <J>T(tk,T)di 

tk-l 


(37) 
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where  x(tk)  is  the  predicted  estimate  vector,  4>  is  the 
state  transition  matrix  der>  'ed  from  the  system  F matrix), 
and  P(tk)  is  the  state  covariance  (prediction  error) 
matrix.  The  superscripts  - and  ♦ denote  before  and  after 
a new  measurement  (update)  is  taken.  The  argument  tk  ^ 
represents  the  sample  time  of  the  preceding  sample.  The 
superscript  T is  the  matrix  transpose  operator.  Since  the 
control  input  is  constant  between  samples,  the  control 
transition  matrix,  T,  is  given  by 

fck 

r(tk'tk-l)  = 1 *(tk,T)L(x)dT  (38) 

fck-l 

where  L is  the  control  distribution  matrix.  The  distur- 
bance distribution  matrix,  G,  describes  which  states  are 
corrupted  by  process  noise  and  the  matrix  Q(t)  is  the 
process  noise  covariance  matrix  given  by  Eq  (7)  where  w(t) 
is  the  process  noise  vector. 

The  Kalman  filter  update  equations  (in  discrete  form) 

are 

x(t£)  = x(tk)  + K(tk)  [z.(tk)  - Hx(t~)  ] (39) 

P(tk)  = P(t”)  - K(tk)H(t")P(t“)  (40) 

and 

T T . 

K(tk)  = P(t")H(tk)  [H(tk)P(t")H(tk)  + R(tk)]"i  (41) 

where  x(tk)  is  the  updated  estimate  (corrected  by  measure- 
ment) , P(t*)  is  the  updated  covariance  (filter  error)  matrix, 
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and  K(t^)  is  the  Kalman  gain  (optimal  weighting)  matrix. 
Variable  :2(t^)  is  the  sampled  measurement  vector,  H is 
the  measurement  matrix  and  the  superscript  -1  above  the 
bracketed  term  in  Eq  (39)  is  the  matrix  inversion  oper- 
ator. The  matrix  R ( ) is  the  measurement  noise  covari- 
ance matrix  given  by  Eq  (8)  where  v(tk)  is  the  measurement 
noise  vector. 

Figure  9 is  a block  diagram  of  a Kalman  filter.  The 
Kalman  filter  model  contains  the  state  transition  matrix 
(derived  from  the  system  F matrix) , the  H matrix,  and  the 
G matrix,  all  of  which  have  been  previously  derived.  In 
addition,  the  noise  covariance  matrices  are  required  and 
are  derived  below. 

Process  Noise  Covariance  Mode 1 . The  covariance  model 
of  the  white  noise  source  driving  the  noise  filter  was 
determined  using  the  variance  of  the  output  of  the  tilt- 
meter  (angular  position)  with  the  platform  in  the  uncon- 
trolled mode.  "The  uncontrolled  mode  is  defined  by  the 
platform  floating  on  the  pneumatic  cylinders  with  no  con- 
trol feedback  from  the  sensors"  (Ref  4:17).  The  uncon- 
trolled mode  is  represented  by  the  system  dynamics  model 
(platform) , the  tiltmeter,  and  the  process  noise  filter. 
The  model  state  matrices,  F^  and  G^,  are  then 

► 
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Figure  9.  Kalman  Filter  (Ref  11:18) 


The  RMS  quiescent  peak-to-peak  amplitude  of  the  con- 
trolled platform  has  been  determind  to  be  ±0.4  arcseconds. 

The  RMS  excursion,  which  can  be  considered  to  represent 
the  standard  deviation  or  one-sigma  value,  is,  therefore, 

0.2  arcseconds.  Since  the  tiltmeter  has  a gain  of  100 
millivolts  per  arcsecond,  the  excursion  corresponds  to 
0.02  volts  output  from  the  tiltmeter.  The  variance  is 
found,  by  squaring  the  one-sigma  value,  to  be  0.0004  volts. 

The  equivalent  process  noise  covariance  model  was 
determined  by  solving  for  Q in  the  steady-state  linear 
covariance  propagation  equation 

P(t)  = 0 = FjP(t)  + P(t)F^  + Gj^OgJ  (43) 

where  P(t)  is  the  time  derivative  of  the  system  covariance. 
The  diagonal  elements  of  P(t)  are  the  variances  of  the 
states.  The  noise  covariance  matrix  is  a one-by-one  matrix 
(scalar)  in  this  model  and  the  p33  element  of  the  covar- 
iance matrix  is  the  variance  of  the  output  of  the  tilt- 
meter (state  x^) , in  volts. 

The  process  noise  covariance  was  found  by  utilizing 
a computer  program  to  integrate  Eq  (43)  and  by  varying  the 
value  of  Q until  the  tiltmeter  variance  approached  the 
value  determined  for  quiescent  excursions.  A value  of 
1.849x10m  for  Q resulted  in  a tiltmeter  variance  of  0.000403. 
This  value  of  Q was  used  in  the  optimal  Kalman  model. 

Measurement  Noise  Covariance  Model.  There  have  been 
no  accurate  models  developed  for  the  noise  characteristics 
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of  the  sensors  employed  in  the  platform  control  system. 

An  estimate,  however,  was  made  by  using  one-half  of  the 

threshold  values  of  the  sensors  as  "rough"  models  of  the 

one-sigma  noise  amplitude.  The  noise  variances  are  the 

one-sigma  values,  squared.  Therefore,  for  the  tiltmeter, 

-4 

the  threshold  is  1.0x10  volts  and  the  variance  is  2.5 

-9  -3 

xlO  volts.  The  threshold  of  the  seismometer  is  1.66x10 

volts  and  the  variance  is  6.889x10  ^ (Ref  12). 

Since  there  are  two  independent  measurement  noise 

sources,  the  measurement  noise  matrix,  R,  is  the  two-by- 

two  matrix  (described  by  Equation  8) . 


R 


2.5x10 


-9 


6.889x10 


-7 


(44) 


Analysis  of  Optimal  Kalman  Filter  Performance 

Kalman  filter  performance  can  be  analyzed,  without 
actually  implementing  the  filter,  by  analyzing  a time 
history  of  the  covariance  of  the  estimates,  P~  and  P+. 

This  is  possible  because  the  covariance  update  and  propa- 
gation equations,  as  well  as  the  associated  Kalman  gain 
equation,  are  not  dependent  on  the  measurement  realiza- 
tions or  the  estimates. 

A program  called  the  "General  Covariance  Analysis 
Program" (GCAP) (Ref  13)  was  used  for  the  filter  performance 
analysis  (and  tuning)  described  in  this  report.  Essentially, 
this  program  generates  a "true"  covariance  time  history 
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for  a filter  at  a given  sampling  rate.  More  specifically, 
the  product  used  most  for  this  study  is  the  one-sigma 
values  (square-roots  of  diagonal  elements)  for  the  system 
state  estimate  errors.  GCAP  is  described  in  more  detail 
in  Appendix  A. 

The  truth  model  represents  the  best  available  model 
of  the  real  world  system.  The  covariance  analysis  of  the 
Kalman  filter  based  on  the  truth  model,  at  a given  sampling 
rate,  represents  a theoretical  performance  bound  at  that  sam- 
pling rate.  This  performance  bound  or  "benchmark"  is  used 
to  evaluate  the  effects  of  varying  the  analysis  parameters, 
i.e.  sampling  frequency,  noise  levels,  sensor  sensitivities, 
and  simplifying  model  reductions. 

The  time  histories  of  the  xl  state  (angular  position, 

8)  and  the  x2  state  (angular  rate,  0)  estimate  errors  (one- 
sigma  values) , for  a sampling  rate  of  200  Hz,  are  presented 
in  Figures  10  and  11,  respectively. 

The  shape  of  the  plots  indicate  that,  after  an  initial 
transient  period,  the  one-sigma  values  for  P+  and  P settle 
to  steady  state  conditions.  This  is  indeeed  the  case  since 
it  cam  be  shown  that  stable  time-invarient  systems  driven 
by  stationary  noises  settle  to  constant  one-sigma  values, 
independent  of  the  initial  state  uncertainties  (Ref  11:11-63) 
For  the  seismic  isolation  platform,  the  noises  are  white 
Gaussian  noise  with  noise  strengths  that  do  not  vary  with 
time  (thus  stationary  statistics)  and  the  system  matrices 
are  time-invariant.  This  result  is  important  since  it 
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Figure  10.  Platform  Tilt  (One-Sigma  Estimate  Error) 
200  Hz  Sample  Rate 
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Figure  11.  Platform  Rate  (One-Sigma  Estimate  Error),  200  Hz.  Sampling  Rate 


permits  the  implementation  of  an  approximate  Kalman  filter 
with  a constant  Kalman  gain 

K = P HTR_1  (45) 

CO 

where  is  the  steady-state  covariance.  The  Kalman  filter 
implementation  is  then  greatly  simplified  (see  Chapter  IV) . 
The  steady  state  one-sigma  values  are  the  values  used  for 
performance  analysis  in  this  study. 

Since  the  optimal  estimates  generated  by  the  Kalman 
filter  are  used  as  inputs  of  the  controllers,  the  design 
criteria  for  acceptable  performance  are  that  the  one-sigma 
errors  in  these  estimates  be  at  least  as  good  as  the  re- 
quired controller  performance  specifications,  i.e.  1.0xl0~3 
arcseconds  for  one-sigma  position  error  and  1.667x10  3 
arcseconds/second  for  one-sigma  rate  errors.  In  addition, 
since  the  largest  error  occurs  for  values  of  P_  (just 
before  measurement  update) , this  error  (prediction  error) 
is  used  as  the  measurement  of  interest. 

The  one-sigma  prediction  errors  for  the  optimal  Kalman 
filter,  with  a 200  Hz  sampling  rate,  are 

9 = 9.05x10*^  arcseconds 

-3  (46) 

9 = 3.37x10  arcseconds/second 

The  position  (tilt)  prediction  error  (9)  is  well  within 
specifications  but  the  rate  prediction  error  is  more  than 
two  orders  of  magnitude  larger  than  required.  Since  this 
Kalman  filter  is  based  on  the  "best"  estimate  of  the 
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physical  system,  the  results  indicate  that  a Kalman  filter 
will  not  provide  the  required  accuracy  with  a 200  Hz 
sampling  rate. 


One  postulated  improvement  is  to  sample  at  a higher 
rate  since,  as  the  sampling  interval  becomes  shorter,  the 
prediction  error  decreases.  '’’o  test  this  postulate,  time 
histories  were  generated  at  various  sampling  rates  from 
the  Nyquist  frequency  (40  Hz)  to  200kHz.  The  one-sigma 
prediction  errors  obtained  from  this  analysis  are  presented 
in  Table  I.  The  plots  of  the  one-sigma  rate  prediction 
errors  are  presented  in  Figures  28  through  36  in  Appendix 
B.  A plot  of  rate  prediction  error  versus  sampling  rate 
is  presented  in  Figure  12. 

With  a three-order-of -magnitude  increase  in  the  sam- 
pling rate  (from  200  Hz  to  200  kHz) , only  slightly  more 
than  a one-order-of-magnitude  decrease  in  the  sampling  rate 
occurs.  Obviously,  a 200  kHz  sampling  is  impossible  to 
implement  and,  in  fact,  a sampling  rate  that  high  would 
tend  to  invalidate  the  white  noise  assumptions.  However, 
by  analyzing  the  filter  at  extreme  sampling  rates,  some 
insight  is  gained  about  performance  bounds  of  the  filter. 
Therefore,  it  is  concluded  that  it  is  not  possible  to  meet 
the  rate  prediction  error  criteria  with  the  platform, 
actuators,  and  sensors  as  presently  configured,  at  any, 
physically  realizable  sampling  rate. 

Since  it  might  be  possible  to  improve  the  prediction 
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Table  I 


Comparison  of  One-Sigma  Prediction 
Errors  at  Different  Sampling  Rates 


Sampling 

Rate 

Angular 

Position  Error 
(arcseconds) 

Angular 

Rate  Error 
(arc seconds/second) 

40 

Hz 

5.7xl0-4 

4 . 5xl0-2 

50 

Hz 

3. IxlO"4 

3 . 3xl0~2 

100 

Hz 

4 . 4xl0-5 

9 . 9x10” 3 

200 

Hz 

9 . 05xl0“6 

3. 37xl0'3 

2 

kHz 

6. 7xl0-7 

5 . 5xl0~4 

20 

kHz 

2.5xl0"7 

2 . 4xl0~4 

200 

kHz 

l.4xl0-7 

9.0xl0~5 

error  by  decreasing  the  error  (noise)  sources,  an  analysis 

of  the  effects  of  decreasing  the  process  noise  or  increasing 
the  sensitivities  of  the  sensors  was  performed.  The 
results  of  that  analysis  (in  the  form  of  one-sigma  pre- 
diction errors)  are  presented  in  Table  II.  The  .plots  of 
the  one-sigma  rate  prediction  error  time  history  are  pre- 
sented in  Figures  28  through  36  in  Appendix  B. 

As  expected,  the  prediction  error  is  most  sensitive 
to  changes  in  the  noise  associated  with  the  angular 
acceleration  measurement  (the  seismometer  sensitivity) , 
but  a one-order-of-magnitude  increase  in  sensitivity  (re- 
duction in  measurement  noise)  does  not  significantly  reduce 


Sampling  Rate  (Hertz) 

Figure  12.  Platform  Rate  Prediction  Error  Versus  Sampling  Rate 


Table  II 
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Effects  of  Decreasing  the  Error 
(noise)  at  200  Hz 

Sources 

Case 

Type  of  One- Sigma 

Improvement  Tilt  Error  (arc- 

seconds) 

One-Sigma  Rate 
Error  (arcseconds/ 
seconds) 

1 

No  Change  9.05x10  ^ 

3. 37xl0-3 

2 

One-order-of-  4.9x10  ^ 

magnitude  re- 
duction of  pro- 
cess noise 

1. 6xl0-3 

3 

One-order-of-  6.5x10  ^ 

magnitude  in- 
crease in  seismom- 
eter sensitivity 

2. 8xl0*3 

4 

One-order-of-  8.4x10  6 

magnitude  in- 
crease in  tilt- 
meter  sensitivity 

3 . Ixl0~3 

5 

Both  Case  3 and  4.1x10  ^ 

Case  4 

1. 9xl0-3 

the  prediction  error.  Even  with  the  sensitivities  of  both 
sensors  improved  by  one  order  of  magnitude,  the  prediction 
error  is  far  above  the  specification.  It  is  concluded, 
then,  that  improvements  in  the  sensors,  in  the  present  con- 
figuration, will  not  bring  the  filter  rate  prediction  error 
down  significantly! 

The  error  is  due,  primarily,  to  the  fact  that  the 
rate  estimate  is  based  on  a position  measurement  and  an 
acceleration  measurement.  With  a direct  measurement  of 
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the  angular  rate,  it  is  possible  that  the  Kalman  filter 
would  provide  estimates  for  the  rate  that  are  accurate 
within  the  specified  levels. 


I 


Since  it  has  been  concluded  that  the  optimal  Kalman 
filter  will  not  meet  the  performance  criteria,  with  the 
system  as  presently  configured,  and,  since  this  study  is 
constrained  to  the  analysis  of  the  system  as  it  presently 
exists,  two  approaches  for  further  analysis  were  considered. 
The  first  approach  is  to  neglect  the  angular  rate  specifi- 
cation and  design  the  simplest  Kalman  filter  that  will 
meet  the  angular  position  specification.  This  is  the 
approach  taken  by  Brunson  (Ref  4) . The  second  approach 
is  to  design  the  Kalman  filter  that  meets  or  exceeds  the 
angular  position  specification  and,  in  addition,  provides 
the  smallest  angular  rate  prediction  error.  It  is  the 
latter  approach  that  was  taken  in  this  study. 

The  analysis  completed  to  this  point  has  not  con- 
sidered the  effects  that  accrue  from  the  implementation 
of  the  Kalman  filter  on  a small,  relatively  slow,  computer 
(PDP-11/03) . The  analysis  has  been  completed  using  a very 
fast,  60-bit  wordlength  machine  and  the  final  implementa- 
tion is  done  on  a slower,  16-kit  wordlength  machine.  The 
filter  performance  is  degraded  further  due  to  the  limita- 
tions of  the  smaller  machine.  These  degradations  are  dis- 
cussed as  part  of  the  algorithm  design  and  implementation 
considerations  (Chapter  IV) . 

As  noted  previously,  the  Kalman  filters  in  this  in- 
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vestigation  can  be  modeled  as  constant  Kalman  gain  filters. 
This  is  a result  of  the  fact  that  the  covariances  (both 
propagated  and  updated)  reach  steady  state  values,  for 
all  of  the  system  states,  after  an  initial  transient  period. 
The  transient  period  is  short,  compared  to  the  time  the 
platform  is  in  use  and,  therefore,  after  an  initial 
"warm-up"  time  for  the  filter  (equal  to  the  transient  time 
for  the  covariances  to  reach  steady  state) , a constant 
Kalman  gain  filter  can  be  used.  The  constant  Kalman  gain 
implementation  greatly  reduces  the  computation  time  re- 
quired by  the  filter  algorithm  because  the  covariance  equa- 
tions and  the  Kalman  gain  equation  (Equations  37,  40  and 
41) , which  are  the  most  time  consuming  computations  in  the 
filter  algorithm,  are  not  computed  as  a part  of  the  real- 
time filter.  The  constant  Kalman  gain  matrix  is  stored, 
in  memory,  in  the  computer  for  use  in  the  estimate  update 
equation  (Equation  39) . In  addition,  a constant  Kalman 
gain  implementation  eliminates  the,  often  severe,  numeric 
difficulties  caused  by  the  covariance  update  equation 
(Equation  40)  that,  normally,  drive  the  wordlength  require- 
ments. The  numeric  difficulties  often  drive  the  implemen- 
tation of  the  Kalman  filter  to  some  type  of  square-root 
form,  e.g.  LJ-D  Factorization,  to  overcome  the  large  word- 
length  requirements.  Since  the  covariance  update  equation 
does  not  drive  the  wordlength  considerations  in  this  in- 
vestigation, an  eigenvalue  test  was  used  to  determine 
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wordlength  requirements. 

In  order  to  determine  what  wordlength  is  required 
to  implement  the  optimal  Kalman  filter,  the  effects  of 
quantization  on  the  system  eigenvalues  were  examined. 

A program  called  STM  (modified  version  of  program  described 
in  Ref  4)  determines  the  normalized  eigenvalue  shifts 
caused  by  quantizing  the  system  to  a finite  wordlength. 

A listing  of  STM  is  presented  in  Appendix  C.  The  eigen- 
values of  the  state  transition  matrix  are  the  roots  of  the 
system  characteristic  equation  (in  discrete  form).  The 
computer  program  STM  computes  the  state  transition  matrix 
for  the  system  and  computes  the  true  eigenvalues  associated 
with  this  state  transition  matrix.  Next,  STM  quantizes 
the  state  transition  matrix  for  various  wordlengths  and 
computes  the  eigenvalues  associated  with  the  quantized 
state  transition  matrices.  The  program  computes  the  dis- 
tance from  the  true  eigenvalues  to  the  unit  circle  and  the 
distance  (shift)  between  the  true  eigenvalues  and  the 
quantized  eigenvalues  for  each  wordlength.  The  criteria 
for  accepting  a given  wordlength  for  implementation  are 
that  the  system  remains  stable  after  quantization  (roots 
inside  unit  circle  on  z-plane)  and  that  the  normalized 
eigenvalue  shift  is  less  than  10  percent  for  each  eigen- 
value at  that  wordlength.  The  normalized  shift,  SQ,  is 
defined  as 
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where  |D  | is  the  magnitude  of  the  distance  (shift)  between 

4 

the  true  eigenvalue  and  the  quantized  eigenvalue  and  J D { 
is  the  magnitude  of  the  distance  from  the  true  eigenvalue 
to  the  unit  circle.  The  wordlength  is  acceptable  only  if 
the  10  percent  criterion  (Ref  4:35)  is  met  by  each  eigen- 
value. In  Figure  13,  the  maximum  normalized  eigenvalue 
shift  is  plotted  for  wordlengths  from  10  bits  to  20  bits. 

The  wordlength  where  all  the  eigenvalue  shifts  are  less 
than  10  percent,  and  thus  the  minimum  acceptable  wordlength 
is  shown  to  be  13  bits.  The  fact  that  the  PDP-11/03  mini- 
computer has  more  bits  than  required  (it  is  a 16-bit  machine) 
allows  more  flexibility  in  scaling  and  lessens  the  effects 
of  overflow  due  to  arithmetic  operations. 

Reduced  Order  Kalman  Filter 

Since  the  optimal  Kalman  filter  is  based  on  a ten-state 
model,  the  number  of  computations  required  to  implement 
the  filter  limits  the  range  of  possible  sampling  rates 
and  it  is  possible  that  a reduced  order  model  (sub-optimal 
Kalman  filter)  might  provide  a smaller  prediction  error 
since  a higher  sampling  rate  could  be  employed. 

In  order  to  investigate  this  possibility,  four  sub- 
optimal  Kalman  filter  models  are  developed  below.  The 
state  matrices  associated  with  the  reduced  models  are  pre- 
sented in  Appendix  D. 

Noise  Filter  Reduction.  The  frequency  response  of 
the  process  noise  filter  is  shown  in  Figure  14.  The  noise 
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filter  provides  a sharp  cutoff  for  frequencies  above 
approximately  20  Hz  (125  radians/second,  see  Equation  24) . 
Since  the  20  Hz  bandwidth  of  the  noise  from  the  process 
noise  filter  (exponentially  correlated)  is  much  greater 
than  the  1.3  Hz  natural  frequency  of  the  platform,  the 
noise  is  assumed  to  be  white  (limited  bandwidth  system 
driven  by  relatively  broadband  noise) . This  simplifying 
assumption  results  in  a reduced  (approximate)  model  com- 
posed of  seven  states.  Since,  as  indicated  by  the  fre- 
quency response,  the  process  noise  filter  attenuated  the 
process  noise  by  approximately  125  dB,  a rough  estimate 
of  the  equivalent  white  noise  for  the  reduced  model  is 

Qeg  = 10~13  QTs  = 1.849X10"1  (48) 

where  Q is  the  estimated  eouivalent  white  noise  strength 
eq 

and  Q,ps  is  the  process  noise  covariance  used  in  the  truth 
model  for  the  system. 

Tiltmeter  Model  Reduction . In  order  to  reduce  the 
tiltmeter  model,  a second-order  approximation  was  developed 
that  has  approximately  the  same  frequency  response  as  the 
tiltmeter  truth  model  (Ref  4:26).  The  form  of  this  approxi- 
mation is 


C _ 4.9(1  - 0.55s)  Volts 

R s2  + 8s  + 49  Arcsecond 


(49) 


A comparison  of  the  frequency  responses  of  the  truth 
model  and  reduced  model  is  presented  in  Figure  15.  As 
indicated  by  the  plot,  the  frequency  response  of  the 
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reduced-order  model  has  the  same  general  shape  as  the 
true  model.  The  responses  appear  to  be  identical  at  very 
low  frequencies  (up  to  approximately  2.5  Hz)  and  differ 
by,  at  most,  1 dB  (beyond  20  Hz) . Figure  16  is  a plot 
of  the  phase  response  of  the  tiltmeter  models.  The 
general  shape  of  the  plots  are  the  same  but  the  phase  re- 
sponse of  the  true  model  lags  the  response  of  reduced- 
order  model.  The  lag  is  a maximum  of  approximately  20 
degrees  (at  2 Hz)  and  is  generally  much  less  than  10  de- 
grees. Since  the  measurement  from  the  tiltmeter  is  con- 
sidered most  accurate  in  the  0-1  Hz  range,  the  reduced- 
order  model  is  an  adequate  representation.  However,  the 
reduced-order  model  will  inject  additional  inaccuracies 
(noise)  into  the  system. 

Seismometer  Model  Reduction.  The  frequency  response 
of  the  seismometer  is  depicted  in  Figure  17.  The  seis- 
mometer exhibits  essentially  no  dynamics  in  the  0-20  Hz 
frequency  band  and  is,  therefore,  approximated  by  a con- 
stant gain  of  one.  This  simplifying  assumption  reduces 
the  system  by  two  states.  Two  reduced  models  were  developed 
using  this  approximation.  First,  in  order  to  investigate 

the  effects  of  the  seismometer  reduction  in  combination 

4 

with  the  true  tiltmeter  model,  a five  state  model  was 
developed.  The  five  state  model  consists  of  the  platform 
dynamics  model  (two  states) , the  true  tiltmeter  model 
(three  states) , the  reduced-order  seismometer  model  (con- 
stant gain  of  1) , and  the  reduced-order  noise  filter  model 
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(constant  strength  white  noise) . Second,  in  order  to  inves- 
tigate the  effects  of  the  combination  of  the  reduced-order 
models,  a four  state  model  was  developed.  The  four  state 
model  consists  of  the  platform  dynamics  (2  states) , the 
reduced-order  tiltmeter  model  (2  states) , the  reduced-order 
seismometer  model  (constant  gain  of  1) , and  the  reduced- 
order  noise  filter  model  (constant  strength  white  noise) . 

In  Chapter  IV,  each  filter  model  (in  combination  with 
the  optimal  controller)  is  analyzed  at  an  appropriate  sam- 
pling rate,  determined  by  the  number  of  calculations 
involved  in  that  filter/controller  implementation. 

Summary 

It  has  been  shown  that,  with  the  seismic  isolation 
platform  as  presently  configured,  optimal  estimation  by 
Kalman  filter  techniques  will  not  reduce  the  uncertainties 
inherent  in  the  system  to  levels  low  enough  to  meet  the 
minimum  accuracy  required  for  successful  optimal  control. 

In  addition,  neither  decreases  in  process  noise,  nor 
increases  in  sensor  sensitivities  significantly  improve 
the  performance  of  the  Kalman  filter.  Possible  significant 
improvement  can  be  achieved  by  augmenting  the  system  with 
a direct  measurement  of  the  angular  rate  of  the  platform. 

Since  the  angular  rate  specification  (1.667x10  5 
arcseconds/second)  cannot  be  met  (the  angular  position 
specification  is  surpassed  by  more  than  two  orders  of 
magnitude) , the  approach  followed  in  this  thesis  is  to 
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design  for  the  best  possible  angular  rate  performance  that 
can  be  expected  for  this  system  configuration.  In  addition 
to  the  optimal  Kalman  filter  model,  four  sub-optimal 
models  were  developed  for  performance  comparisons.  The 
approach  is  to  select  and  implement  the  filter  model  that 
provides  the  best  rate  prediction  performance  (smallest 
error) . The  selected  Kalman  filter  is  combined  with  the 
optimal  controller  developed  in  Chapter  III. 


i 

' 

1 
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III.  Controller  Development 

The  development  of  the  controller  segment  is  based 
on  discrete  models  of  the  platform  and  the  actuators 
(pneumatic  cylinders  and  shakers) . A control  segment 
model  was  proposed  by  Burkhart  (Ref  1) . The  model  was 
developed  for  a sampling  rate  of  200  Hz  and  tested  with 
a 1.25  foot-pound  step  input.  This  chapter  contains  a 
brief  discussion  of  the  models  and  techniques  employed 
in  the  development  of  the  controller  segment.  In  addi- 
tion, pertinent  results  stated  by  Burkhart  are  included 
as  background  information.  In  Chapter  IV,  a control  seg- 
ment is  developed,  in  general  form,  based  on  the  sampling 
rate  dictated  by  the  analysis  of  the  various  filter/con- 
troller combinations. 

Separation  of  Controller  Segment 

The  design  of  the  controller  segment  was  divided  into 
two  tasks  (Fig  8) . Each  control  task  is  associated  with 
one  of  the  two  types  of  actuators  used  in  controlling  the 
angular  motion  about  the  horizontal  axis.  The  justifica- 
tion for  the  separation  of  the  controller  segment  is  based, 
primarily,  on  the  dynamic  characteristics  of  the  actuators. 
The  pneumatic  actuators  are  slow  (time  constant  of  20 
seconds)  and  are  actually  part  of  the  platform  support 
system.  Due  to  the  pneumatic  actuators'  dynamic  response, 
they  are  employed  in  a position  feedback  loop  to  counter  any 
torque  imbalances  resulting  from  unsymmetrical  loading 
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of  the  platform.  The  shakers  have  a much  faster  dynamic 
response  (natural  frequency  is  26  Hz)  and  are  employed 
in  an  optimal  state-feedback  control  (optimal  regulator) 
loop,  to  regulate  the  effects  of  the  environmental  dis- 
turbances. In  both  cases,  the  control  loops  were  developed 
with  the  assumption  that  they  are  receiving  perfect  infor- 
mation about  the  system  states  from  the  Kalman  filter,  i.e. 
exact  knowledge  of  the  entire  state. 

The  optimal  regulator  is  designed  to  regulate  the 
pneumatic  loop.  The  states  associated  with  the  pneumatic 
actuator  are  the  contents  of  registers  in  the  computer 
(derived  from  the  pneumatic  loop  compensator  algorithm) 
and  are  known  exactly.  These  states  are  incorporated  into 
the  optimal  control  problem  through  the  use  of  an  augmented 
state  space  model  that  is  used  in  deriving  the  control  law. 
Assuming  the  exact  values  of  the  states  (from  the  Kalman 
filter  and  from  the  pneumatic  loop) , the  deterministic 
discrete-time  optimal  controller  will  minimize  the  discrete 
performance  index  given  in  Eq  (3) (Ref  15:502).  This  is, 
again,  based  on  the  "forced  separation"  concept  described 
in  Chapter  I. 

Continuous  System  Models 

The  development  of  the  controllers  is  based  on  the 
models  of  those  components  used  in  controlling  angular 
motion  about  the  horizontal  axis.  The  sensor  models  are 
not  included  in  the  controller  segment  development  since 
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they  are  part  of  the  estimator  problem  ("forced  separation") 
The  models  required  include  the  seismic  isolation  platform 
model,  the  pneumatic  actuator  model,  the  shaker  actuator 
model,  and  a model  of  the  environmental  disturbances  (pro- 
cess noise) . 

Seismic  Isolation  Platform.  The  transfer  function  for 
the  seismic  isolation  platform  was  presented  in  Eqs  (7)  and 
(8)  (Chapter  II) . 

Pneumatic  Actuator.  The  pneumatic  actuator  transfer 
function  is  (Ref  1:17) 


Hpls) 


s + .05 


Foot-pound 

Volt 


Shaker  Actuator . The  shaker  actuator  transfer  function 


is  (Ref  1:18) 


k (i) 

H (S)  = 

s + 2^susS 


Foot-pounds 
2 Amps 


where 


«s  = °*7 

wg  = 157  radians/second 

and  kg  = 70  foot-pounds/amp 

Environmental  Disturbance  (Noise)  Mode 1 . To  evaluate 
the  performance  of  the  controller,  a deterministic  model 
of  the  environmental  disturbance  (called  process  noise  in 
Chapter  II)  was  developed.  A torque  step  function,  ut(t), 
directly  applied  to  the  top  surface  of  the  concrete  block 
was  selected  since  it  is  readily  modeled  in  both  the  z-domain 
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and  the  s-domain.  There  is  some  question  as  to  the  weight 
of  the  step-input  that  more  closely  represents  the  pro- 
cess noise.  The  peak  overshoot  of  the  time  response  of 
the  platform  model,  to  a step  input,  approximates  the  aver- 
age excursion  of  the  uncontrolled  platform  due  to  the  envi- 
ronmental disturbances.  Burkhart's  model  is  based  on  an 
average  excursion  of  0.1  arcseconds , derived  from  Reference 
12:2.  The  average  excursion  used  to  determine  the  process 
noise  model  for  the  Kalman  filter  was  0.2  arcseconds  (Ref 
8) . Since  no  estimate  of  the  quiescent  platform  oscillatory 
excursion  appears  to  be  more  acceptable  than  the  others, 
the  "worst  case"  model  (0.2  arcseconds)  is  used  in  this 
investigation.  However,  for  the  purpose  of  describing  the 
results  of  Burkhart’s  work,  the  0.1  arcsecond  excursion 
model  is  employed. 

Using  the  computer  program  TOTAL,  the  weight  of  the 
representative  step  input  was  determined  to  be  1.25  foot- 
pounds. As  seen  in  Figure  18,  a 1.25  foot-pound  step  input 
produces  a peak  overshoot  of  0.102  arcseconds. 

Figure  18  illustrates  the  relationships  of  the  torques 
and  transfer  function  associated  with  the  actuators  and 
the  process  noise. 

Discrete  System  Models 

Two  discrete  models  were  derived  from  the  continuous 
system  transfer  functions.  A z-domain  function,  0(z),  of 
angular  position  of  the  platform,  with  a step  input,  was 
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ISOLATION  PLATFORM  TIME  RESPONSE  TO  1.25  STEP 


Figure  18.  Isolation  Platform  Time  Response  to  1.25  Foot-Pound 

Step  Input 


Figure  19.  Block  Diagram  of  Platform  and  Actuators 

(Ref  1:21) 

derived  for  use  in  designing  the  controller  for  the  pneu- 
matic control  loop.  A discrete  state  representation  of  the 
platform  and  the  actuators  was  developed  for  use  in  the 
designing  of  the  optimal  controller  in  the  regulator  loop. 

Position  Feedback  Loop.  A block  diagram  of  the  pneu- 
matic loop  is  shown  in  Figure  20. 

The  pneumatic  loop  includes,  along  with  the  platform 
(G^fs))  and  the  pneumatic  actuator  (H  (s)  ) , an  impulse 
sampler,  a computer  algorithm  D(z)  and  a digital-to-analog 
converter  (DAC) . The  impulse  sampler  in  the  model  prepresents 
the  estimation  process  completed  by  the  Kalman  filter  (de- 
noted 9*(s)).  The  DAC  is  represented  by  a zero-order  hold 
(ZOH) , H#(s).  The  ZOH  holds  the  output  of  the  DAC  constant 
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Pneumatic  Zero-  Digital 

Actuator  Order-Hold  Controller 


1 


I 


Figure  20.  Block  Diagram  of  Pneumatic  Loop  (Ref  1:22) 

between  samples  and  has  the  transfer  function 


Hc(s) 


(52) 


where  T is  the  sampling  period. 

The  computer  algorithm,  D(z),  is  a digital  computer 
program  (digital  controller  or  compensator)  that  receives 
the  estimated  angular  position  from  the  Kalman  filter  and 
produces  a control  signal,  C^(s),  that  is,  in  turn,  sent 
to  the  DAC  and  converted  to  analog  form  for  input  to  the 
pneumatic  actuator. 

The  location  of  the  computer  in  the  feedback  loop 
prohibits  the  reduction  of  the  block  diagram  to  the  z-domain 
transfer  function  0(2'/ut(z).  Instead,  the  approach  is  to 
solve  for  an  expression  of  0(z).  First 


. 
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E(s)  = ut(s)  - Hp(s)Ho(s)CD(s) 


0 (s)  = ut(s)Gb(s)  - Gb(s)Hp(s)Ho(s)cJ(s) 


And  taking  the  z-transform 


0(2)  = utGb(z)  - GbHpHo(z)CD(z) 


Since 


CD(z)  = D(z)0(z) 


Then 

u G.(z) 

6(2)  = 1 + GbHpHQ(z)D(z)  (57) 

The  z-transforms  u.G,  (z)  and  G,  H H (z)  must  be  deter- 

t d d p o 

mined.  Since  the  design  objective  is  to  meet  the  specified 
response  characteristics  for  a step  input  ut(s)Gb(s)  is 
evaluated  with 


ut(s)  = 1/s 


where 


Bz  + Cz 

Z [u.  (s)  G.  (s)  ] = 5— 

D (z-1)  (z2+£z+e) 


Ae  aTsin(bT  - 0) 


= -2e~aTcos (bT) 


-2aT 


kb(B  - a + 1) 
kb(a  + e) 


I 


The  coefficients,  evaluated  at  T = 5 msec,  are 
B = 2.6915839xl0-5 

C = 2 . 68844  53xl0-5 

6 = 1.9952834 

e = 0.9965061 2 

The  z-transform  for  G.  (s)H  (s)H  (s) , since  H (s)  is 

p o o 

a zero-order  hold,  is 


G.H  H ( z ) 
b p o 


z-1 

z' 


(s)hp 

s 


(65) 


Now 


(s) Hp (s) 
s 


4 3 2 

Jz  + Lz  + Mz  + Nz 

(z-1) (z-6) (z2+6z+e) 


(66) 


where 


J = 

G + k,  - E 
b 

(67) 

L = 

8kb  - BE  + 

E - 6k.  + H - 6G  - G 
b 

(68) 

M = 

ekb  - eE  + 

6 Bk.  + EB  + 6G  - H6  + H 
b 

(69) 

N = 

6ekb  + eE  + 

6G 

(70) 

The  coefficients,  calculated  at  T = 5 msec  are 
J = 2. 1439340 xlO-17 

L = 2. 2435920xl0-9 

M = 8. 9654109xl0-9 

N = 2. 2393891xl0“9 

6 = 0.99975003 

To  verify  the  discrete  models,  an  initial  value  veri- 
fication was  employed.  For  verification,  the  initial  values 
of  the  continuous  time  functions  were  compared  to  the  initial 
values  for  the  corresponding  discrete-time  functions.  All 
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the  functions  had  initial  values  of  zero  except  the  dis- 
crete function  G^M^z)  which  had  an  initial  value  equal 
to  the  value  of  the  coefficient  J.  Since  J is  several 
orders-of-magnitude  smaller  than  the  other  coefficients, 
the  term  involving  J,  i.e. , the  z4  term,  was  deleted. 

This  deletion  is  justified  since,  with  single-precision 
16-bit  wordlength  computations,  J is  too  small  to  be  repre- 
sented and  would,  therefore,  be  rounded  to  zero.  The  re- 


sulting function  is 

G,  (s) H (s) 

Z [ — —E ] 


Lz3  + MZ2  + Nz 


(z-1) (z-6) (z  — Bz+e) 


To  obtain  G.H  H (z) , Equation  71  is  substituted  in  Equation 
— P o 

65  to  yield 


G.H  H (z)  = 

b p o 


Lz  + Mz  + N 
(z-6)  (z2-Bz+e) 


Finally,  substituting  Equations  72  and  59  into  Equation  57 


yields 


9(2)  = 


Bz  + Cz 
(z-1) (z2+Bz+6) 


Lz  +MZ+N 

(z-6) (z2+6z+e) 


Discrete  State  Space  Model.  Since  the  zero-order  hold 
device  maintains  the  control  inputs  constant  over  the  sampling 
period,  the  continuous  state  space  model  is  transformed  to 


a discrete  state  space  model  by  the  following  transforma- 


tion : 


$ = e 
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r “ /e(T'T)BdT  (75) 


o 

where  $ is  the  discrete  state  transition  matrix,  F is  the 
continuous  plant  matrix,  T is  the  sampling  interval,  T is 
the  control  transition  matrix  (discrete  control  matrix) , and 
B is  the  continuous  control  matrix.  The  continuous  state 
equation  is 

x(t)  = Fx(t)  + Lu(t)  (76) 

with  u(t)  for  te[kT,  (k+l)T]  and  the  discrete  state  equa- 
tion is 

x (k  +1)  = $x(k)  + ru(k)  (77) 

Using  phase  variable  form  state  equations  to  represent 
the  individual  transfer  functions,  the  resulting  state 
matrices  are 


F 


0 


1 0 


0 0 


-49  -.7 


1 1725430  0 


0 0 -.05  0 

0 0 0 0 


0 

1 


0 


0 -24649  -219.8 


(78) 


and 


0 


0 


0 0 


0 


L 


.05  0 


0 0 


0 


1 


(79)  I 
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where 


0 = 2.156x1(t) 

0 = 2.156x2(t) 


and 


u(t)  = 


up(t) 

us(t) 

ut(t) 


The  general  form  of  the  transition  matrices  is 


(80) 

(81) 


(82) 


*11 

*12 

*13 

*14 

*15 

*21 

*22 

*23 

*24 

*25 

0 

0 

*33 

0 

0 

0 

0 

*43 

*44 

*45 

0 

0 

*53 

*54 

*55 

ru 

r!2 

ri3 

r21 

r22 

r23 

= 

r31 

0 

0 

0 

r 

0 

0 

r 

0 

mm 

(83) 


(84) 


The  numerical  values  of  the  elements  of  the  matrices 
is  determined  by  the  relationships  described  in  Equations 
74  and  75  and  is  a function  of  the  sampling  interval,  T. 


Pneumatic  Controller  (Compensator) , D(z) 

The  criteria  used  in  designing  the  pneumatic  loop  corn- 
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pensator,  D(z),  are  that  the  system  remains  stable  with 
D(z)  in  the  feedback  loop,  that  D(z)  drives  the  platform 
to  zero  steady-state  angular  position,  and  that  the  tran- 
sient response  (peak  overshoot  and  settling  time) , to  a 
1.25  foot-pound  input,  of  the  controlled  platform  is  better 
than  that  of  the  uncontrolled  platform  (depicted  in  Figure 


The  steady  state  response  of  6(z)  was  analyzed  by 


letting 


D ( z ) = 


ND(z) 

d^T 


in  Equation  73  and  employing  the  final  value  theorem 


(t-*-°°)  = lim(z-l)  [ 
z-1 


z (Az+B) (z-6) 0D (z) 
(z-1) P ( z / 


where 


P (z)  = DD(z)  [z3  + (0  + <5)z2  + (e  -6  0)  z -<5e] 

+ (Lz2  + Mz  + N)Nd(z)  (87) 

To  insure  that  the  steady  state  value  of  9(z)  goes 
to  zero,  D ( z ) must  have  a factor  (z-1)  in  the  denominator 
(Ref  16:290) . 

Root  locus  techniques  were  used  to  examine  the  stability 
and  transient  response  of  various  compensators  having  a 
factor  of  (z-1)  in  their  denominators.  Using  the  root 
loci  of  the  expression,  D ( z) G^HpHQ (t)  (Equation  58),  for 
various  D(z)'s  the  characteristic  values  of  the  pneumatic 
loop  were  determined.  By  varying  the  gain  and  the  location 
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and  number  of  poles  and  zeros  in  the  compensator,  a form 
for  D ( z ) was  chosen  that  produces  a stable  system  and  a 
fast  transient  response.  The  selected  form  for  the  dis- 
crete controller  for  the  pneumatic  loop  was 

, . 7500  (z-Q.  9 9900)  (z 2-2 . 0150z+l . 0158 ) 

(z-l. 0)  (z2-0 .9752)  (88) 

The  response  of  the  system,  with  the  compensator,  was 
determined  by  evaluating  the  difference  equation  (trans- 
formed from  9(z),  Equation  73)  with  the  above  D(z)  in  the 
feedback  loop  (Ref  1:52).  As  depicted  in  Figure  21,  the 
steady-state  response  for  the  controlled  system  is  zero 
(versus  0.055  arcseconds  for  the  uncontrolled  platform, 
Figure  17)  and  the  peak  overshoot  is  slightly  less  than 
that  for  the  uncontrolled  platform  (0.097  arcseconds  versus 
0.102  arcseconds).  However,  the  settling  time  is  slightly 
greater  for  the  compensated  system  (approximately  11 
seconds  versus  10.86  seconds).  The  D(z)  described  above 
was  designed  based  on  a sampling  rate  of  200  Hz.  It  is 
assumed,  for  analysis  purposes,  that  the  general  form  of 
D (z) (Equation  88)  is  maintained  through  the  change  in 
sampling  rate.  This  assumption  is  based  on  the  fact  that 
changes  in  the  sampling  rate  will  cause  changes  in  the  lo- 
cations of  the  roots  of  the  discrete  system  models  (plat- 
form and  pneumatic  actuator) , but  will  not  change  the 
order  of  these  models.  Therefore,  the  general  form  of  the 
compensator  that  is  used  in  Chapter  IV  for  determining  the 
approximate  sampling  rate  is 
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Figure  21.  Graph  of  Difference  Equation  for  9 (kT) 

(Ref  1:53) 

2 

k (z-a) (z  -bz+c) 

D(z)  = — 5 (89) 

(z-1. 0) ( z -dz ) 

Burkhart  determined  that,  with  the  form  of  D(z)  in 
Equation  88,  the  time  responses  (to  a 1.25  ft-lb  step-input) 
for  the  discrete  compensator  and  for  a thirteen-bit  digitized 
version  of  the  compensator  agreed  within  three  significant 
digits  after  25  seconds.  It  was  concluded  that  the  thirteen 
bit  controller  shows  no  significant  degradation  of  perfor- 
mance (Ref  1:57). 

Regulator  Loop . The  regulator  (shaker)  loop  is  used 
to  provide  further  damping  of  the  pneumatic  loop  and  also 
to  provide  control  of  the  angular  rate  of  the  platform. 

The  shaker  receives  its  control  input  from  a state  feed- 
back control  law  (optimal  controller)  that  is  based  on  an 
augmented  system  state  model  consisting  of  the  discrete  state 
space  model  of  the  inertial  instrument  test  platform  and 


actuators  (Equations  83  and  84)  combined  with  a discrete 
state  space  model  of  the  pneumatic  compensator. 


A state  space  model  of  the  pneumatic  actuator  is  de 


rived  by  factoring  Equation  88  into 


rw  x _C(z)  _ 7500  z-0.999  z2-2 . 015z+l. 0158 

D[Z)  TTzT  z (2-1.0)  (z-0.975) 


and  forming  a state  space  representation  for  each  term 
using  partial  fraction  expansion.  After  combining  the 
state  representations  for  each  term,  the  state  space  modi 
of  D(z)  is  (see  Figure  19) 


O 

o 

o 

l 

■ — 

1 

xD(k+l)  = 

-0.999  1.0  0 

^(k)  + 

1 

-0.999  0 0.975 

_ — 

1 

C (k)  = 7500{  [-0.999  3. 36324xl0-2  -7.36324x10 

+ 0 (JO  } 


( 


] xD(k) 
( 


After  augmenting  the  system  model  (Equations  83  and 
84)  with  the  compensator  model  (Equation  91  and  92) , the 
discrete  pneumatic  loop  state  space  model  becomes  (in  ge 
eral  form) 
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where  $ is  the  pneumatic  loop  state  transition  matrix 
and  Tp  is  the  pneumatic  loop  control  transition  matrix. 


The  objective  of  the  deterministic  discrete  optimal 
controller  is  described  by  Equations  3 and  4 in  Chapter 

A 

I except  that,  due  to  the  forced  separation,  x(tk)  (the 
state  estimate  from  the  Kalman  filter)  is  replaced  by 
x(t^)  for  design  purposes.  Thus,  the  design  of  the  con- 
trol law,  C(t^)  is  assumed  independent  of  the  stochastic 
properties  of  the  system. 

The  quadratic  performance  index,  J (in  Equation  4) , 
can  be  interpreted  as  "system  error  plus  control  effort" 
measure  of  performance  that  uses  a tradeoff  between  system 
error,  represented  by  quadratic  term  involving  the  V 
matrix,  and  control  effort,  represented  by  the  quadratic 
control  "intensity"  term  involving  the  U matrix  (Ref  17: 
326)  . 

The  procedure  used  to  solve  for  C(t^)  was  adapted  from 
Linear  Optimal  Control  Systems  (Ref  15:502)  and  consists 


' 


of  solving  the  equations 

c(tk)  = {u  + Tp  [v  + p (k+i)  ] r }-1rT  [v  + p(k+i)]$  05) 

and 

P (k)  = $T[V  + P (k+1)  ] [4>  - FC(k)]  (96) 

backward  in  time  from  the  terminal  condition 

P(n)  = Vf  (97) 

It  has  been  demonstrated  that  the  solution  of  Equa- 
tions 95  through  97  results  in  a steady-state  solution  for 
C(tk)  that  is  independent  of  the  terminal  condition.  The 
resulting  control  law  is  time  invariant  and  asymptotically 
stable  (Ref  1:62). 

The  solution  of  the  optimal  control  law  does  not 
guarantee  that  the  design  specifications  will  be  met  (Ref 
1:61).  The  optimal  Gontrol  law  must  be  solved  by  the 
iterative  process  of  selecting  various  values  for  the 
weighting  matrices. V,  V^,  and  U,  and  evaluating  the  re- 
sulting control  law  to  determine  if  it  meets  the  design 
specifications.  The  resulting  control  law  vector,  in  gen- 
eral form,  is 

T 

(98) 


'll 

12 

13 

14 

15 

16 

17 

18 
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The  results  of  Burkharts  analysis  (Fig  21)  indicate 
that,  for  the  control  law  based  on  a sampling  rate  of 
200  Hz,  the  angular  rate  specification  is  not  met,  al- 
though the  system  does  settle  to  within  the  design  speci- 
fication (1.667xl0~5  arcseconds/second)  within  0.08  seconds. 
It  was  concluded  that  the  angular  rate  specification  could 
not  be  met,  at  any  physically  realizable  sampling  rate, 
by  the  designed  optimal  controller,  due  to  the  fact  that 
the  angular  rate  of  the  inertial  test  platform  was  already 
three  orders-of-magnitude  greater  than  the  design  specifi- 
cation at  the  end  of  the  first  sample  period  (0.05  seconds) 

(Ref  1:68).  The  optimal  controller  controlled  the  angular 

-4 

position  to  less  than  2.9x10  arcseconds  (much  better  than 
the  design  specification  of  1x10  3 arcseconds) . 

Summary 

In  this  chapter,  the  general  forms  of  the  control  equa- 
tions are  developed.  The  pneumatic  loop  controller  uses 
angular  position  feedback  to  generate  the  control  signal 
to  the  pneumatic  actuator.  The  matrix  form  of  the  con- 
trollers are  converted  to  scalar  equation  form  to  alleviate 
the  inefficiency  of  matrix  operations  with  sparse  matrices 
and  also  to  allow  more  flexibility  in  the  ordering  of  the 
computations.  The  resulting  pneumatic  controller  scalar 
equations  are  (from  Equation  91  and  92) . 

x6(k+l)  = 2.156  xx(k)  (99) 
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(100) 


up  = CD(k)  = kDtuP(k)  + x6(k+l)] 
x?(k+l)  = -0.999  xfi (k)  + x?(k)  + Xg(k+1)  (101) 

xg(k+l)  = -0.999  x& (k)  + 0.975  Xg  (k)  +xg  (k+1)  (102) 

Up (k+1)  = 0.999  Xg  (k+1)  + 3. 36xl0"2x? (k+1) 

- 7 . 36x10 _2Xg (k+1)  (103) 

The  optimal  regulator  uses  state  feedback  to  derive 
a control  law  that  produces  the  control  signals  to  the 
shaker  actuator.  The  scalar  equations  representing  the 
optimal  regulator  algorithm  are  (from  Equations  3 and  98) 

ug(k)  = c11x1(k)  + c21x2(k)  + c31x3(k) 

+ c42x4 (k)  + c5ix5 (k)  + Ug(k)(104) 

us(k+l)  = c61Xg(k+l)  + c71x7(k+l)  + Cg^Xg (k+1)  (105) 

In  Chapter  IV,  these  scalar  equations  are  used  to 
determine  a rough  estimate  of  the  computation  times  in- 
volved in  the  implementation  of  various  filter/controller 
pairs.  Although  the  results  of  Burkharts  investigation 
indicate  that  the  angular  rate  specification  cannot  be  met, 
even  when  the  input  from  the  cascaded  Kalman  filter  is 
assumed  to  be  within  design  specifications,  the  angular 
rate  specification  is  not  neglected  in  the  remainder  of 
this  investigation.  As  mentioned  in  Chapter  II,  the  design 
approach  is  to  provide  improved  angular  rate  control  and 
to  meet  or  exceed  the  angular  position  design  specification. 


76 


Based  on  the  results  of  Burkhart's  investigation,  the 
optimal  control  approach  warrants  further  investigation 
as  a method  of  isolating  the  inertial  instrument  test 
platform  at  FJSRL. 


77 


IV.  Selection  of  Kalman  Filter  Model 


All  of  the  model  development  in  the  previous  chapters 
was  based  on  a sampling  rate  of  200  Hz.  The  selection 
of  this  sampling  rate  was  based  on  engineering  judgement 
i.e. , it  is  five  times  the  Nyquist  frequency.  In  reality, 
since  the  computations  involved  in  executing  the  optimal 
estimation  and  control  algorithm  take  a finite  length  of 
time,  the  sampling  rate  will  be  limited  by  the  instruction 
execution  (arithmetic,  store,  shift,  etc.)  times  of  the 
PDP-11/03  Minicomputer.  In  this  chapter,  the  Kalman  filter 
models  developed  in  Chapter  II  are  combined,  in  turn,  with 
the  general  model  of  the  controllers  developed  in  Chapter 
III  to  determine  the  approximate  maximum  sampling  rate  for 
each  combination.  Each  filter  is  then  tuned,  using  the 
General  Covariance  Analysis  Program  (GCAP) , to  determine 
an  expected  performance  bound  for  that  filter.  The  tuned 
filter  performances  are  compared  and  a "best  cut"  filter 
model  is  selected  as  a baseline  for  future  implementation. 

In  addition,  some  implementation  considerations  are  pre- 
sented as  background  for  possible  "follow-on"  investigation. 

General  Form  of  the  Algorithm 

The  implementation  of  the  optimal  estimation  and  con- 
trol algorithm,  developed  in  this  investigation,  is  sep- 
arated into  five  tasks: 

- Input 


- Output 


- Supervisor 

- Pneumatic  Control  Loop 

- Optimal  Regulator  Loop 


\ 


| 

J 


“ Kalman  Filter 

Input.  The  input  task  is  to  move  the  current  measure- 
ment from  the  Analog-to-Digital  Converter  (ADC)  to  a 
storage  register  within  the  computer.  Assuming  the  use 
of  the  PDP-11/03  compatable  ADC  (Model  ADV-ll-A) , the 
measurement  enters  the  computer,  after  conversion,  in  off- 
set binary  form.  The  input  routine  checks  for  overflow 
or  underflow  and  the  measurement  data  is  converted  to 
two's  complement  form  since  the  arithmetic  operations  used 
by  this  algorithm  are  computed  using  two's  complement 
arithmetic. 

Output.  The  output  routine  is  essentially  the  mirror 
image  of  the  input  routine.  After  all  computations  on  con- 
trol data  are  complete,  the  output  routine  converts  the 
data  from  two's  complement  form  to  offset  binary  form, 
checks  for  overflow  or  underflow,  and  moves  the  data  from 
a storage  register  to  the  Digital-to-Analog  Converter  (DAC) . 

Supervisor.  The  supervisor  routine  is  used  to  con- 
trol the  sequencing  of  the  other  tasks.  In  addition, 
various  "housekeeping"  tasks  (i.e.,  scaling,  error  checking, 
etc)  can  be  accomplished  in  this  routine. 

The  controllers  and  the  Kalman  filter  are  implemented 
in  the  remaining  subroutines  and  are  discussed  in  more 
detail  below. 
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Figure  22.  Implementation  Tasks  for  Optimal  Estimation 
and  Control  Algorithm. 


Figure  22  is  a block  diagram  representation  of  the 
tasks  involved  in  the  optimal  estimation  and  control  al- 
gorithm. 


Computation  Times 


Based  on  the  general  form  of  the  algorithm  described 
above,  the  total  computation  time  for  the  optimal  estimation 


and  control  algorithm  can  be  described  as 


( (Tj  + Tq  + Tp  + Tr  + Tr)*1.1) 


(106) 


where 


= total  computation  time 
= input  computation  time 


Tq  = output  computation  time 

Tp  = pneumatic  control  loop  computation  time 

Tr  = optimal  regulator  computation  time 

T = Kalman  filter  computation  time 
K 

The  1.1  multiplier  represents  a 10%  overhead  that  is  added 
to  represent  the  supervisory  operations. 

Since  there  are  actually  two  optimal  estimation  and 
control  systems  involved  in  isolating  the  platform  (two 
horizontal  axes  through  the  center  of  the  platform) , the 
total  approximate  computation  time,  Tar,  will  be 

Tac  = 2TC  (107) 

For  purposes  of  evaluating  the  sampling  rates  for 
the  various  Kalman  filter  models,  the  computation  times 
for  the  input,  output,  pneumatic  controller,  and  optimal 
regulator  are  considered  to  be  the  same  for  each  filter 
case. 

The  instruction  execution  times,  for  the  PDP-11/03 
Minicomputer,  for  the  instructions  employed  in  the  optimal 
estimation  and  control  algorithm  are  presented  in  Table 
III  (Ref  18) 

All  arithmetic  operations  are  assumed  to  be  single- 
precision fixed-point  operations.  Floating-point  and/or 
double  precision  arithmetic  operations  would  provide  more 
accuracy  for  representing  the  terms,  but  it  is  felt  that 
the  increases  in  computation  times  involved  more  than  off- 
set the  benefits  derived.  All  coefficients  are  normalized 
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Table  III 


PDP-11/03  Instruction  Execution  Times 


Assembly  Execution  Time 

Instruction  Language  Code  (microseconds) 


Move 

MOV 

3.5 

Compare 

CMP 

6.65 

Fixed-Point 

Subtraction 

SUB 

7.7 

Fixed-Point 

Addition 

ADD 

7.7 

Fixed-Point 

Multiply 

MUL 

64.0 

Arithmetic  : 

Shift  Combined 

ASHC 

15.0 

(by  scaling)  so  that  they  are  represented  as  fractions. 

In  this  way,  the  multiplication  operation  will  always 
result  in  a fraction  and,  therefore,  multiplication  over- 
flows will  be  avoided.  To  illustrate  the  result,  the 
following  simple  example  is  presented 


•L2 

• 510 

X •12 

= 

•510 

(108) 

.oi2 

• 2510 

where  . 1 2 is  the  binary  representation  of  .5  (decimal). 
The  PDP-11/03  multiplication  operation  will  result  in  the 
following  product 


(109) 
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which  is  obviously  an  incorrect  result.  In  order  to  cor- 
rect the  result,  an  Arithmetic  Shift  Combined  (AHSC)  in- 
struction is  used  to  shift  the  result  right  by  one  position 
(with  a zero  fill) , resulting  in  a correct  representation 
of  the  product.  An  AHSC  instruction  is  provided  with  each 
multiplication  instruction  in  this  investigation.  As  in- 
dicated in  Table  III,  this  will  increase  the  multiplication 
computation  time  by  15  microseconds  for  each  operation, 
resulting  in  a total  multiplication  time  of  79  microseconds. 

Input/Output  Computation  Times . The  input  routine  con- 
sists of  one  move  instruction,  two  compare  instructions 
(checking  for  overflow/ underflow) , and  one  fixed-point  sub- 
traction instruction  (to  convert  from  offset-binary  to  two's 
complement  form)  and  results  in  a total  computation  time  of 
approximately  24.5  microseconds.  The  output  routine  is 
essentially  the  mirror  image  of  the  input  routine  (except 
one  fixed-point  addition  instead  of  subtraction  to  convert 
from  two's  complement  to  offset  binary  form)-  and  also  re- 
sults in  a total  computation  time  of  approximately  24.5 
microseconds.  The  analog-to-digital  and  digital-to-analog 
conversion  times  are  not  included  in  the  calculation  of  the 
I/O  computation  times  since  it  is  assumed  that  an  interrupt 
scheme  will  be  used  for  controlling  the  ADC,  thus  permitting 
simultaneous  conversion  and  filter/controller  algorithm 
execution,  and  it  is  also  assumed  that  the  DAC  can  complete 
its  conversion  process  without  interferring  (based  on  Model 


ADV-11— A Converter)  with  filter/controller  algorithm  ex- 
ecution . 

Pneumatic  Control  Loop  Algorithm.  The  pneumatic  con- 
trol loop  algorithm  is  represented  by  Equations  99  through 
103  in  Chapter  III..  The  algorithm  consists  of  a series  of 
multiplications  and  additions.  It  is  assumed  that  the  common 
(constant)  gain  operation,  KD,  in  Equation  100,  is  accom- 
plished using  an  external  analog  amplifer.  An  example 
of  the  steps  involved  in  the  execution  of  the  algorithm  is 
(using  Equation  101) 


MOV 

TERMl , R2 

(Move  -0.999  into  Register  2) 

MUL 

x6  , R2 

(Multiply  x, (k)  by  -0.999  and 
store  the  product  in  Registers 

2 and  3) 

MOV 

TERM2 , R4 

(Move  0.975  into  Register  4) 

MULT 

Xg  ,R4 

(Multiply  Xg (k)  by  0.975  and 
store  the  product  in  Register  4 

ADD 

R4  ,R2 

(Add  the  products  and  store 
result  in  Register  2) 

ADD 

Xg (k+1) ,R2 

(Add  Xg  (k+1)  to  previous  sum, 
store  result  in  Register  2) 

MOV 

R2 , Xg (k+1) 

(Move  the  result  to  memory 
location  of  Xg(k+1) 

The  total  pneumatic  control  algorithm  requires  13  move 
instructions,  6 fixed-point  additions,  and  5 fixed-point 
multiplications,  resulting  in  a total  computation  time  of 
approximately  565.4  microseconds. 

Optimal  Regulator  Algorithm.  The  optimal  regulator 
algorithm  is  represented  by  Equations  104  and  105.  In 
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[ 

a manner  similar  to  the  precedure  used  for  the  pneumatic 
control  algorithm,  the  computation  time  for  the  optimal 
regulator  loop  was  determined  to  be  approximately  694.4 
microseconds . 

The  total  computation  time  for  the  input,  output,  and 
controllers  is  approximately  1308.8  microseconds.  This 
computation  time  is  assumed  constant  for  all  filter/con- 
troller combinations. 

Kalman  Filter  Algorithms . The  Kalman  filter  algorithms 
are  represented  by  Equations  34  and  38.  In  order  to  allow 
for  control  inputs,  each  Kalman  filter  model  is  augmented 
with  the  state  associated  with  the  pneumatic  actuator.  The 
augmented  state  represents  the  output  torque  generated  by 
the  pneumatic  actuators.  Although  the  augmented  models  are 
not  used  during  the  covariance  analysis  to  determine  filter 
performance  ("forced  separation"),  they  are  included  in  the 
implemented  algorithm  and  do,  therefore,  affect  the  com- 
putation time  for  each  filter.  Therefore,  the  state  matrices 
referred  to  in  the  Kalman  filter  computation  time  deter- 
mination are  the  state  matrices  associated  with  the  aug- 
mented filter  models.  These  models  are  presented  in  Appendix 
E.  To  determine  the  computation  times  for  the  estimate 
propagation  equation  (Equation  34)  for  each  filter,  the 
general  form  of  the  state  transition  matrix  (4>)  and  the 
discrete  control  matrix  (T)  was  derived  (using  the  STM 
program  with  a sampling  rate  of  200  Hz)  for  each  filter 
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model.  By  converting  the  matrix  equations  to  scalar  form, 
and  applying  the  procedure  employed  for  determing  the  con- 
troller computation  times,  the  approximate  computation 
time  for  the  propagation  equation  was  determined  for  each 
filter. 

To  determine  the  computation  time  for  the  estimate 
update  equation  (Equation  38) , the  matrix  equation  was 
again  converted  to  scalar  form  which  resulted  in  the  set 
of  equations  represented  below. 

A , A 

x^tt^)  = x^t”)  + kil  [Residual  1]  + ki2  [Residual  2]  (110) 

The  residual  terms  represent  the  difference  between 
the  current  measurement,  z(t^),  and  the  predicted  measure- 

A 

ment,  Hx(t£),  and  are  common  to  each  scalar,  state  estimate 
update  equation  for  a given  filter  model,  therefore,  the 
residuals  are  computed  once  for  each  sample  period.  As 
mentioned  previously,  a constant  Kalman  gain  implementation 
is  employed,  and  therefore,  the  covariance  equations  and 
the  Kalman  gain  equation  are  not  computed  on-line.  The 
computation  times  for  the  estimate  update  equation,  for 
each  filter  model,  was  determined  in  a manner  similar  to 
that  described  for  the  propagation  equation. 

The  total  approximate  computation  times,  TAC,  and 
the  corresponding  maximum  sampling  rates  for  the  various 
filter/controller  combinations  are  presented  in  Table  IV. 
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Table  IV 


Filter/Controller  Maximum  Sampling  Rates 


Filter 

T (Eqs  106 
and  107) 
(milliseconds) 

Maximum 
Sampling  Rate 

Optimal  (10  state) 

22.17 

45.1 

Seven-state 

15.63 

64.0 

Six-state 

12.04 

83.0 

Five-state 

10.95 

91.4 

Four-state 

9.46 

105.7 

Kalman  Filter  Performance  Evaluations 

Using  the  sampling  rates  indicated  in  Table  III,  the 
General  Covariance  Analysis  Program  was  employed  to  deter- 
mine the  expected  performance  bound  for  each  Kalman  filter 
model.  The  object  of  this  evaluation  is  to  investigate 
the  possibility  that  one  of  the  suboptimal  Kalman  filters, 
based  on  the  increased  sample  rate  (relative  to  the  optimal 
filter)  permitted  by  the  reduction  in  the  number  of  com- 
putations required  by  that  reduced-order  model,  provides 
better  performance  than  the  optimal  Kalman  filter  sampled 
at  the  rate  dictated  by  the  large  number  of  computations 
involved.  The  performance  measures  used  as  criteria  for 
this  evaluation  are  the  platform  rate  prediction  errors 
which  are  produced  as  an  output  from  the  GCAP  computer 
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program.  Again,  the  criteria  are  that  the  angular  posi- 
tion specification  is  met  or  exceeded  and  that  the  angular 
rate  prediction  error  is  as  small  as  possible. 

Optimal  Ka lman  Filter.  Figures  23  and  24  are  plots 
of  the  time  histories  of  the  one-sigma  errors  for  the  plat- 
form angular  position  and  angular  rate,  respectively,  from 

the  optimal  Kalman  filter  sampled  at  45.1  Hz.  The  steady- 
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state  prediction  errors  are  2.53x10  Arcseconds  for  angular 
position  and  2.92xl0-2  Arcseconds/Second  for  angular  rate. 

Seven-State  Kalman  Filter.  A seven-state  suboptimal 
Kalman  filter  model  was  developed  in  Chapter  II  based  on  a 
simplifying  assumption  which  approximated  the  third-order 
noise  filter  model  with  a roughly  equivalent  white  noise 
model.  This  approximation  inserts  more  uncertainty  into 
the  filter  model.  The  Kalman  filter  tuning  process  consists 
of  inserting  pseudo-noise  into  the  reduced-order  model, 
by  increasing  the  strengths  of  the  process  noise  (matrix  Q) 
and/or  the  measurement  noises  (matrix  R) , until  the  true 
system  error  approximates  the  error  generated  by  the  re- 
duced-order filter.  By  varying  the  white  noise  strength 
represented  by  the  Q matrix,  the  equivalent  white  noise 
driving  the  platform,  that  tuned  the  filter,  was  deter- 
mined to  be  approximately  3.25.  This  is  not  a unique  solu- 
tion, since  other  values  for  Q and  R (in  combination)  will 
tune  the  filter.  The  resulting  time  histories  for  the 
angular  position  and  angular  rate  are  presented  in  Figures 
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25  and  26.  The  steady-state  prediction  errors  are  3.34x10 
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Figure  24.  One-sigma  Optimal  Filter  Rate  Error  (45.1  Hz) 


Figure  25.  One-Sigma  Seven-State  Filter  Tilt  Error  (64  Hz) 


PLAT 


_ 2 

Arcseconds  for  angular  position  and  3.78x10  Arcseconds/ 

Second  for  angular  rate. 

Other  Filter  Models.  In  the  process  of  tuning  the 
six-state,  the  five- state , and  the  four-state  suboptimal 
Kalman  filters,  the  preliminary  results  indicated  that  the 
error  performances  from  those  filters  were  far  worse  (on 
the  order  of  200  Arcseconds/Second  for  the  one-sigma  rate 
prediction  error)  than  the  performance  indicated  for  the 
higher-order  filters.  Based  on  these  preliminary  results, 
it  was  felt  that  the  time-consuming  task  of  tuning  these 
reduced-order  filters  would  be  little,  if  any,  benefit  to 

I 1 

the  performance  evaluation. 

Filter  Selection.  The  performance  evaluation  indicated 
that,  based  on  the  maximum  permissible  sampling  rates,  the 
optimal  Kalman  filter  provides  better  performance  tha.i  any 
of  the  suboptimal  filters.  However,  the  performance  of  the 
seven-state  filter  approaches  that  of  the  optimal  filter 
and  should  not  be  eliminated  from  consideration.  The  in- 
creased sampling  rate  possible  for  the  seven-state  filter 
will  tend  to  improve  the  performance  of  the  associated  con- 
trollers, and  it  is  possible  that  the  improvement  in  con- 
troller performance  due  to  the  increased  sampling  rate, 
for  the  suboptimal  filter,  is  greater  than  the  improve- 
ment due  to  better  Kalman  filter  performance  from  the 
optimal  filter.  In  order  to  investigate  this  possibility, 
two  pneumatic  control  loop  compensators  and  two  associated 
optimal  regulators  would  be  designed,  one  controller  pair 
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based  on  a sampling  rate  of  45.1  Hz  and  the  other  based 
on  a sampling  rate  of  64  Hz.  Each  controller  pair  would 
be  combined,  in  turn,  with  the  associated  Kalman  filter 
and  the  total  performance  of  the  optimal  estimation  and 
control  system  would  be  evaluated.  Due  to  time  constraints, 
this  performance  analysis  was  not  undertaken  in  this  in- 
vestigation but  is  recommended  as  part  of  a follow-on  in- 
vestigation . 

A mulitarized  version  of  the  LSI-11  hardware  multiplier 
is  available  that  reduces  the  multiplication  execution  time 
from  64  microseconds  to  12.2  microseconds.  If  the  faster 
multiplier  is  used,  the  maximum  sampling  rates  for  the 
various  filter  model  implementations  are  as  indicated  in 
Table  IV. 

As  expected,  since  the  multiplication  operation  is  the 
most  time  consuming  operation,  the  large  reduction  is 
multiplication  execution  time  results  in  a significant 
increase  in  maximum  sampling  rates.  Since  the  increase 
in  sampling  rates  is  significant,  a re-evaluation  of  the 
filters  at  these  sampling  rates  is  warranted  and  is  recom- 
mended as  part  of  a follow-on  investigation. 

Implementation  Considerations 

Of  prime  importance  in  implementing  the  optimal  esti- 
mation and  control  algorithm  is  the  time  lag  between  the 
time  the  sampled  measurements  are  performed  and  the  time 
the  associated  control  signal  is  output.  In  order  to 
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Table  V 


Maximum  Sampling  Rates  Using  High-Speed  Multiplier 


Filter 

Sampling  Rate  (Hz) 

Optimal 

111 

Seven-state 

147 

Six-state 

218 

Five-state 

237 

Four-state 

239 

reduce  the  magnitude  of  the  state  propagation  (prediction) 
error  generated  by  the  state  propagation  equation  (Equa- 
tion 34) , the  computation  lag  time  should  be  minimal.  To 
ensure  the  least  possible  lag  time,  only  those  computations 
that  actually  require  the  updated  measurement  and  those 
directly  associated  with  the  control  output  should  be  com- 
puted between  measurement  update  and  control  output.  A 
proposed  sequence  of  computations  is  presented  in  Figure 
27. 

The  two  optimal  estimation  and  control  systems  associ- 
ated with  the  horizontal  axes  through  the  center  of  gravity 
are  assumed  to  be  independent  control  systems.  This 
assumption  is  based  on  the  assumption  of  decoupled  modes 
of  motion  mentioned  in  Chapter  I.  Since  each  of  these 


Sampled  Measure-  Sampled  Measure- 
ments at  kT  -<ents  at  (k+l/2)T 

| Control  system  l| IX 


2 3 4 


12  3 4 


Sampled 
Measurements 
at  1 (k+l)T 

Time 


1.  Sample  and  Input  Measurement 

2.  Update  state  estimate 

3.  Compute  and  output  pneumatic  control  signal 

4.  Compute  and  output  optimal  regulator  control 

5.  Start  background  computations 


After  the  sampled  measurements  are  input,  the  next 
step  is  to  determine  the  updated  state  estimates  using 


Equation  38.  Using  the  forced  separation  concept,  since 
some  of  the  states  estimated  by  the  Kalman  filter  are  not 
used  in  the  control  output  computations,  only  those  states 
that  are  used  in  the  control  algorithm  need  to  be  updated 
at  this  time,  the  remaining  state  updates  will  be  computed 
as  background  computations,  immediately  before  the  state 
propagation  computations. 

The  next  step  is  to  compute  and  output  the  new  control 
signals,  beginning  with  the  pneumatic  control  signal.  In 
the  pneumatic  compensator  development  discussed  in  this 
thesis,  the  state  variables  were  chosen  so  that  the  state 
associated  with  the  input  to  the  pneumatic  loop  compensator 
(Xg)  is  proportional  to  the  platform  angular  position  (x^) 
as  represented  in  Equation  80.  By  choosing  the  state  vari- 
ables so  that  the  state  xg  is  equal  to  state  x^,  the  multi- 
plication required  in  Equation  99  would  be  eliminated  and 
the  pneumatic  control  signal  could  be  computed  simply  by 
solving  Equation  100.  In  addition,  by  using  an  external 
analog  amplifier  to  perform  the  multiplication  associated 
with  the  term  KD  (after  the  control  signal  is  output) , the 
computation  lag  time  is  reduced. 

Finally,  using  Equation  104,  the  optimal  regulator 
control  signal  is  generated  and  output. 

The  remaining  equations  are  solved  as  background  com- 
putations in  the  time  remaining  until  the  samples  are  taken 
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tions,  the  optimal  Kalman  filter  provides  the  best  per- 
formance among  the  filter  models  developed.  It  was  deter- 


mined that  the  optimal  Kalman  filter  can  be  implemented 

with  a maximum  sampling  rate  of  approximately  45  Hz  and, 

at  this  sampling  rate,  the  angular  position  specification 
-4 

is  met  (2.53x10  Arcseconds)  but  the  angular  rate  specifi- 
. -2 

cation  is  not  met  (2.92x10  Arcseconds/Second) . The 
performance  of  the  filter  is  improved  if  the  sampling  rate 
is  increased,  and  various  methods  of  increasing  the  sampling 
rate  were  discussed.  In  addition,  a proposed  sequence 
for  algorithm  computation  was  presented  for  consideration 
for  follow-on  investigations. 
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V.  Conclusions  and  Recommendations 

The  purpose  of  this  investigation  was  to  determine 
the  expected  performance  derived  from  implementing  linear 
stochastic  estimation  and  control  techniques  to  the  problem 
of  actively  controlling  the  inertial  instrument  test  plat- 
form at  FJSRL.  The  study  was  divided  into  six  areas: 

1.  The  development  of  a state  space  system  model 
from  which  an  optimal  Kalman  filter  was  developed  and 
evaluated  using  a covariance  analysis  tehcnique. 

2.  The  development  of  four  suboptimal  (reduced-order) 
Kalman  filter  models. 

3.  The  development  of  general  models  for  a pneumatic 
(position-feedback)  control  loop  compensator  and  an  optimal 
(state-feedback)  regulator. 

4.  The  determination  of  the  approximate  maximum  sampling 
rates  for  each  Kalman  filter  model,  in  turn,  in  combination 
with  the  general  models  of  the  controller  algorithms. 

5.  The  tuning  and  performance  evaluation  of  each  Kalman 
filter  model,  in  turn,  based  on  the  appropriate  maximum 
sampling  rate  for  that  model. 

6 . A proposed  sequence  for  the  execution  of  the 
optimal  estimation  and  control  algorithm  computations. 

Conclusions 

The  specific  conclusions  derived  from  this  investi- 
gation are: 

1.  As  presently  configured,  the  platform  angular  rate 
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uncertainty  specification  (1.667x10  5 Arcseconds/Second) 
cannot  be  met  at  any  physically  realizable  sampling  rate. 

The  failure  to  meet  this  specification  is  due,  primarily, 
to  the  fact  there  is  no  direct  measurement  of  platform 
angular  rate. 

-3 

2.  The  platform  angular  position  specification  (1x10 

Arcseconds)  can  be  met  by  the  optimal  Kalman  filter  (2.53 
>4 

xlO  Arcseconds)  with  a sampling  rate  of  45  Hz.  In  addi- 
tion, the  angular  position  specification  can  be  met  (Ref  1: 
75)  based  on  a sampling  rate  of  200  Hz  and  the  assumption 
that  the  Kalman  filter  provides  estimations  that  are  accu- 
rate within  the  angular  positon  specification.  Further 
investigation  is  required  to  determine  if  the  angular 
position  specification  can  be  met  by  the  optimal  Kalman 
filter  in  cascade  with  controllers  based  on  sampling  rate 
of  45  Hz. 

3.  A suboptimal  (seven-state)  Kalman  filter  model, 
sampled  at  64  Hz,  meets  the  angular  position  specification 

(3.34xl0~4  Arcseconds)  but  fails  to  meet  the  angular  rate 

_o 

specification  (3.78x10  Arcseconds/Second).  Further  in- 
vestigation into  the  performance  of  this  filter  cascaded 
with  the  appropriate  controllers  is  warranted  due  to  the 
possible  benefits  to  be  derived  by  the  increased  sampling 
rate. 

Recommendations 

Five  recommendations  have  resulted  from  this  study: 


100 


1.  A covariance  analysis  of  the  performance  of  the 
system,  with  the  addition  of  an  angular  rate  sensor  model, 
should  be  completed  to  determine  if  the  improvement  in 
performance,  due  to  the  rate  sensor,  warrants  the  inclu- 
sion of  such  a sensor. 

2.  The  stochastic  noise  models  used  in  this  investi- 
gation were  primarily  based  on  engineering  judgement  and 
not  on  any  realistic  empirical  models.  Further  study  should 
be  directed  toward  the  development  of  process  noise  and 
sensor  error  models. 

3.  Further  study  should  be  directed  toward  the  develop- 
ment of  accurate  models  of  all  the  system  components, 
especially  the  inertial  test  platform  (center  of  gravity, 
homogeneity  of  the  structure,  resonances,  bending  modes, 
etc. ) . 

4.  In  order  to  increase  the  maximum  sampling  rates, 
two  recommendations  are  proposed: 

By  employing  the  militarized  version  of  the  LSI-11 
hardware  multiplier  the  multiplication  computation  time 
is  decreased  from  64  microseconds  to  approximately  12 
microseconds  and  substantial  increases  in  sampling  rates 
are  derived. 

- Even  more  substantial  improvement  can  be  derived 
by  utilizing  several  microprogrammable  microprocessors 
to  perform  the  computations.  With  this  approach,  it  is 
possible  to  accomplish  several  operations  simultaneously 
(in  parallel) , resulting  in  increased  sampling  rates. 
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5.  Further  testing  and  evaluation  of  the  optimal 


l 


estimator  and  the  controllers  (in  cascade) , based  on 
physically  realizable  sampling  rates,  should  be  accom- 
plished using  digitally  implemented  (on  the  PDP-11/03 
minicomputer)  algorithms  and  simulated  noise  disturbances. 
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Appendix  A 


General  Covariance  Analysis  Program 

Performance  analyses  were  performed  using  covariance 
analysis  techniques  implemented  on  the  General  Covariance 
Analysis  Program  (GCAP) . Since  the  covariance  equations 
and  the  Kalman  gain  equation  (equations  35,  39,  and  40) 
can  be  computed  independent  of  any  measurements,  it  is 
possible  to  perform  a performance  analysis  of  a Kalman 
filter  design  without  actually  simulating  the  sampling 
process.  The  GCAP  employs  the  covariance  equations  and 
the  Kalman  gain  equation  to  generate  the  statistics  (in 
the  form  of  one-sigma  time  histories)  of  a filter  design, 
directly. 

Two  mathematic  system  models  are  used  by  the  GCAP. 

One  model,  called  the  "truth  model",  represents  the  best 
available  system  model  (used  to  develop  the  optimal  Kalman 
filter)  and  the  other  model,  called  the  filter  model,  is 
a reduced-order  model  used  to  develop  a suboptimal  Kalman 
filter.  As  the  order  of  the  model  is  decreased,  modeling 
errors  are  introduced  into  the  design  and  the  Kalman  filter 
based  upon  such  a reduced-order  model  will  provide  poorer 
performance  than  the  optimal  filter. 

In  GCAP,  both  truth  model  and  filter  model  covariance 
matrices  are  propagated  in  time  using  a fourth-order  Runge- 
Kutta  numerical  integration  formula.  At  each  sample  measure- 
ment update  time,  the  filter  covariance  update  is  performed 


(Equation  39)  and  an  optimal  Kalman  gain  is  computed  based 
on  the  filter  model.  The  Kalman  gain  equation  becomes 


W - VV^fVVf  * v1 


(111) 


w 

P(tk) 

H. 


where 

Optimal  Kalman  filter  (based  on  filter  model) 

Filter  covariance  prediction  matrix 

._  — ' Filter  measurement  matrix 

F 

Rp  * Filter  measurement  noise  matrix 
the  superscript  T is  the  matrix  transpose  operator, 
and  the  superscript  -1  is  the  matrix  inversion  operator. 

The  truth  model  covariance  update  is  computed  using 
the  Kalman  gain  matrix  derived  from  the  filter  model.  The 
magnitude  of  the  truth  model  covariance  elements  (diagonal 
elements)  represent  the  "true  system"  one-sigma  errors. 

The  equation  describing  the  system  update  is 


W - HS,Vtk)(I-MKF<VHS) 


* «VVRsW 

where 

I ■ Identity  matrix 
Hg  ■ System  measurement  matrix 
Rs  - System  measurement  noise  matrix 
and  M ■ is  the  transformation  matrix 

M * I-  - 

0 


(112) 


(113) 


used  to  augment  the  K_  matrix  with  null  elements  in  order 

r 
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to  make  the  matrices  conform  for  arithmetic  matrix  oper- 
ations . 

The  tuning  process  involves  varying  the  tuning  param- 
eters (Qp  and  Rp)  until  the  covariance  performance  of  the 
truth  model  is  at  least  a good  as  the  covariance  performance 
of  the  filter  model.  Varying  the  tuning  parameters  refers 
to  the  process  of  adding  "pseudo-noise"  to  the  filter 
model  to  compensate  for  the  modeling  errors  produced  by 
reducing  the  amount  of  information  known  by  the  system 
(reducing  the  number  of  states) . 


Appendix  B 


The  following  plots  are  one-sigma  error  time  his- 
tories, generated  by  the  General  Covariance  Analysis 
Program,  resulting  from  the  sensitivity  analysis  of 
the  optimal  Kalman  filter,  performed  in  Chapter  II. 
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Figure  29.  Optimal  Filter  Rate  Prediction  Error:  Sampling  Rate  is  50  Hz. 


Figure  31.  Optimal  Filter  Rate  Prediction  Error:  Sampling  Rate  is  2 kHz 


Figure  32.  Optimal  Filter  Rate  Prediction  Error:  Sampling  Rate  is  20  kHz. 
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Figure  33.  Optimal  Filter  (200  Hz)  Low  Process  Noise 
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Figure  35.  Optimal  Filter  (200  Hz)  Increased  Seismometer  Sensitivity. 


Figure  36.  Optimal  Filter  (200  Hz)  Increased  Tiltmeter  and  Seismometer  Sensitivity. 


Appendix  C 


Program  STM  - Quantized  Eigenvalue  Analyses  Program 


A listing  of  the  computer  program  developed  (modified 
from  Ref  3)  to  analyze  finite  wordlength  effects,  based 
on  eigeavalue  shift  criteria,  is  presented  in  this  appendix. 
The  control  cards  required  to  execute  the  program 

are 


ATTACH , A , M3  7 1SUB , 10-AFIT. 

ATTACH , BIMSL , 10-LIBRARY , SN-ASD . 
ATTACH , AFIT , AF I TSUB ROUTINES , ID-AFIT . 
FTN. 

LDSET (LIB-A/B/AFIT) . 

LGO. 


The  input  record  file  consists  of  the  data  card  de- 
scribed in  the  program  and  input  cards  containing  the  state 
matrices  F and  L (in  that  order)  by  rows.  The  format  for 
the  data  cards  is  shown  below  (example  is  for  2x9  matrix) 


Column  10 

Card  1 - Dimension  2 

Cards  2 to  5 - Data  0 

(Right  Justified)  0 

-49 

7 


20 

9 

1 

2.56 

0 

8 


30 

Gi 

1 


40 

0 

2 


50 

0 

4 


60  70 

0 2.156 

5 6 
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FIND  SCALE  FACTOR  THAT  MAKES  ALL  STATE  TRANSITION  MATRIX  ELEMENTS 
LESS  THAN  ONE 
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PRINT  10'*, N 


PRINT*, “ QUANTIZED  MATRIX 
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COMPUTE  AND  PRINT  OISCFETE  CONTROL  MATRIX  (CONTROL  TRANSITION 
MATRIX),  B. 
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Appendix  D 


Sub-Optimal  Filter  Models 
Seven-State  Model 
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